show code
knitr::opts_chunk$set(echo = TRUE)knitr::opts_chunk$set(echo = TRUE)Data and R codes available: github, then it will be submited to zenodo
Inclusion of species diversity of recruiting plants: item I.C, with rarefaction and Proportion of Biotic/Abiotic Seed Dispersal and Planted/Unplanted species.
Table with statistical details from the models (file "itatinga_NUC_ModelSummaries.xlsx")
Fixed LiDAR figure (SE of the mean)
Fixed plots with empty recruits
Considering only small recruits in species richness rarefaction curves
Analyses separating Baccharis + Vernonanthura
path <- getwd()
version _
platform x86_64-w64-mingw32
arch x86_64
os mingw32
crt ucrt
system x86_64, mingw32
status
major 4
minor 3.3
year 2024
month 02
day 29
svn rev 86002
language R
version.string R version 4.3.3 (2024-02-29 ucrt)
nickname Angel Food Cake
if(!require(pacman)){install.packages("pacman")}Carregando pacotes exigidos: pacman
pkg_list <-c("tidyverse", "lme4", "DHARMa", "emmeans", "multcomp", "iNEXT",
"car", "performance", "cowplot", "MuMIn", "ggpattern")
pacman::p_load(pkg_list, character.only = TRUE)trees_clean: data of trees (planted/not-planted)
trees_clean <- readxl::read_xlsx(paste0(path, "/Itatinga_Nuc_trees.xlsx"))
trees_clean <- trees_clean |>
dplyr::mutate(scientificName = dplyr::if_else(is.na(scientificName), "Unclassified", scientificName))grass_cover: grass cover per plot (%)
grass_cover <- readxl::read_xlsx(paste0(path, "/Itatinga_Nuc_grasscover.xlsx"))lidar: Canopy Height Metric (CHM_mean, CHM_SD, and CHM_CV), Canopy cover (Canopy_cover), Rugosity (Rugosity_SD, and Rugosity_CV), Leaf Area Index (LAI_mean, LAI_SD, and LAI_CV), and LAIunder (LAIunder_mean, LAIunder_SD, LAIunder_CV).
lidar2022_wide <- readxl::read_xlsx(
paste0(path, "/tab_metrics_LiDAR_2022_v2.xlsx")) |>
dplyr::select(Trat:LAIunder_SD) |>
tidyr::separate(Trat, into = c("Bloco", "Perc_cobertura","Tratamento"), sep = "_") |>
dplyr::mutate(Tratamento = dplyr::if_else(is.na(Tratamento), "controle", Tratamento),
Tratamento = paste(Perc_cobertura, Tratamento, sep="_"),
ano = 2022) |>
dplyr::select(-Perc_cobertura) New names:
• `` -> `...1`
Warning: Expected 3 pieces. Missing pieces filled with `NA` in 5 rows [5, 10, 15, 20,
25].
lidar2024_wide <- readxl::read_xlsx(
paste0(path, "/tab_metrics_LiDAR_2024_v2.xlsx"))|>
dplyr::select(Trat:LAIunder_SD) |>
tidyr::separate(Trat, into = c("Bloco", "Perc_cobertura","Tratamento"), sep = "_") |>
dplyr::mutate(Tratamento = dplyr::if_else(is.na(Tratamento), "controle", Tratamento),
Tratamento = paste(Perc_cobertura, Tratamento, sep="_"),
ano = 2024) |>
dplyr::select(-Perc_cobertura)New names:
• `` -> `...1`
Warning: Expected 3 pieces. Missing pieces filled with `NA` in 5 rows [5, 10, 15, 20,
25].
lidar_wide <- dplyr::bind_rows(lidar2022_wide, lidar2024_wide) |> # 40 x 12
dplyr::mutate(
Tratamento = factor(Tratamento,
levels = c("25_tabuleiro", "25_faixa", "50_tabuleiro", "50_faixa", "100_controle"),
labels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%")))two tables:
trees_clean |> dplyr::select(RP, family, scientificName) |>
unique()|>
dplyr::mutate(valor=TRUE)|>
tidyr::pivot_wider(names_from = RP, values_from = valor, values_fill = FALSE) |>
dplyr::arrange(desc(planted), family, scientificName) # |> writexl::write_xlsx("tab1.xlsx")# A tibble: 108 × 4
family scientificName planted regenerant
<chr> <chr> <lgl> <lgl>
1 Anacardiaceae Astronium graveolens TRUE FALSE
2 Annonaceae Annona cacans TRUE TRUE
3 Bignoniaceae Handroanthus impetiginosus TRUE FALSE
4 Bignoniaceae Handroanthus ochraceus TRUE TRUE
5 Boraginaceae Cordia sellowiana TRUE TRUE
6 Euphorbiaceae Croton floribundus TRUE TRUE
7 Euphorbiaceae Croton urucurana TRUE FALSE
8 Euphorbiaceae Mabea fistulifera TRUE TRUE
9 Fabaceae Anadenanthera colubrina var. cebil TRUE TRUE
10 Fabaceae Apuleia leiocarpa TRUE FALSE
# ℹ 98 more rows
trees_clean |> dplyr::select(RP, family, scientificName) |>
unique() |> with(table(RP))RP
planted regenerant
37 86
trees_clean |> with(table(RP, scientificName)) |>
as.data.frame() |>
dplyr::filter(Freq>0) |> # writexl::write_xlsx("tab_spp_plantio.xlsx")
dplyr::filter(RP=="planted") RP scientificName Freq
1 planted Anadenanthera colubrina var. cebil 204
2 planted Annona cacans 4
3 planted Apuleia leiocarpa 11
4 planted Astronium graveolens 5
5 planted Cariniana estrellensis 2
6 planted Cecropia pachystachya 117
7 planted Cedrela fissilis 28
8 planted Ceiba speciosa 54
9 planted Citharexylum myrianthum 16
10 planted Copaifera langsdorffii 24
11 planted Cordia sellowiana 10
12 planted Croton floribundus 831
13 planted Croton urucurana 726
14 planted Enterolobium contortisiliquum 77
15 planted Erythrina mulungu 32
16 planted Eugenia uniflora 9
17 planted Ficus guaranitica 82
18 planted Gallesia integrifolia 50
19 planted Genipa americana 3
20 planted Guazuma ulmifolia 323
21 planted Handroanthus impetiginosus 30
22 planted Handroanthus ochraceus 20
23 planted Hymenaea courbaril 12
24 planted Inga edulis 85
25 planted Inga vera 741
26 planted Lafoensia pacari 34
27 planted Luehea divaricata 48
28 planted Mabea fistulifera 20
29 planted Machaerium acutifolium 17
30 planted Myroxylon peruiferum 2
31 planted Peltophorum dubium 29
32 planted Platypodium elegans 7
33 planted Psidium myrtoides 2
34 planted Pterogyne nitens 3
35 planted Senegalia polyphylla 474
36 planted Senna macranthera 1
37 planted Unclassified 372
Reading scientific names
survival_planted2024 <- trees_clean |> # 5416 stems
dplyr::filter(ano==2024)|>
dplyr::select(bloco:arv, Hm, scientificName, RP) |> # 3,199 stems
unique() |> # 1487
dplyr::filter(RP=="planted") |> # 1135
dplyr::group_by(tratamento,bloco, parcela)|>
dplyr::summarise(n = dplyr::n(), .groups = "drop")|>
dplyr::arrange(desc(n))24 seedlings per plot
survival_planted2024 |>
with(table(bloco, parcela, by=tratamento)), , by = 100
parcela
bloco 1 2 3 4
1 1 1 1 1
2 1 1 1 1
3 1 1 1 1
4 1 1 1 1
5 1 1 1 1
, , by = F25
parcela
bloco 1 2 3 4
1 1 1 1 1
2 1 1 1 1
3 1 1 1 1
4 1 1 1 1
5 1 1 1 1
, , by = F50
parcela
bloco 1 2 3 4
1 1 1 1 1
2 1 1 1 1
3 1 1 1 1
4 1 1 1 1
5 1 1 1 1
, , by = N25
parcela
bloco 1 2 3 4
1 1 1 1 1
2 1 1 1 1
3 1 1 1 1
4 1 1 1 1
5 1 1 1 1
, , by = N50
parcela
bloco 1 2 3 4
1 1 1 1 1
2 1 1 1 1
3 1 1 1 1
4 1 1 1 1
5 1 1 1 1
survival_planted2024 |> dplyr::arrange(tratamento, bloco, parcela)# A tibble: 100 × 4
tratamento bloco parcela n
<chr> <chr> <chr> <int>
1 100 1 1 17
2 100 1 2 15
3 100 1 3 12
4 100 1 4 19
5 100 2 1 11
6 100 2 2 11
7 100 2 3 17
8 100 2 4 13
9 100 3 1 13
10 100 3 2 16
# ℹ 90 more rows
mod_surv1 <- lme4::glmer(cbind(n, 24 - n) ~ tratamento + (1|bloco), data = survival_planted2024, family = binomial(link = "logit"))
mod_surv_null <- lme4::glmer(cbind(n, 24 - n) ~ 0 +(1|bloco), data=survival_planted2024, family = binomial(link = "logit"))
AIC(mod_surv1, mod_surv_null) df AIC
mod_surv1 6 554.2907
mod_surv_null 1 564.5393
par(mfrow=c(1,2))
DHARMa::plotResiduals(mod_surv1)
DHARMa::plotQQunif(mod_surv1)Marginal \(R^2\) / Conditional \(R^2\)
MuMIn::r.squaredGLMM(mod_surv1)Warning: the null model is only correct if all the variables it uses are identical
to those used in fitting the original model.
R2m R2c
theoretical 0.1291593 0.4869551
delta 0.1160836 0.4376569
0.129/ 0.487[1] 0.2648871
Which it is 0.265
sjPlot::tab_model(mod_surv1)| cbind(n,24-n) | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 1.10 | 0.79 – 1.52 | 0.580 |
| tratamento [F25] | 0.80 | 0.62 – 1.03 | 0.090 |
| tratamento [F50] | 0.91 | 0.70 – 1.18 | 0.472 |
| tratamento [N25] | 0.58 | 0.45 – 0.75 | <0.001 |
| tratamento [N50] | 0.84 | 0.65 – 1.09 | 0.192 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 bloco | 0.10 | ||
| ICC | 0.03 | ||
| N bloco | 5 | ||
| Observations | 100 | ||
| Marginal R2 / Conditional R2 | 0.010 / 0.038 | ||
broom.mixed::tidy(mod_surv1) |> dplyr::mutate(dplyr::across(4:7, round, 3)) Registered S3 method overwritten by 'broom':
method from
nobs.multinom MuMIn
Warning: There was 1 warning in `dplyr::mutate()`.
ℹ In argument: `dplyr::across(4:7, round, 3)`.
Caused by warning:
! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
Supply arguments directly to `.fns` through an anonymous function instead.
# Previously
across(a:b, mean, na.rm = TRUE)
# Now
across(a:b, \(x) mean(x, na.rm = TRUE))
# A tibble: 6 × 7
effect group term estimate std.error statistic p.value
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 0.092 0.166 0.553 0.58
2 fixed <NA> tratamentoF25 -0.222 0.131 -1.70 0.09
3 fixed <NA> tratamentoF50 -0.094 0.131 -0.719 0.472
4 fixed <NA> tratamentoN25 -0.544 0.132 -4.12 0
5 fixed <NA> tratamentoN50 -0.171 0.131 -1.31 0.192
6 ran_pars bloco sd__(Intercept) 0.309 NA NA NA
summary(mod_surv1)Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: cbind(n, 24 - n) ~ tratamento + (1 | bloco)
Data: survival_planted2024
AIC BIC logLik -2*log(L) df.resid
554.3 569.9 -271.1 542.3 94
Scaled residuals:
Min 1Q Median 3Q Max
-3.0363 -0.8992 -0.1404 1.0640 3.6310
Random effects:
Groups Name Variance Std.Dev.
bloco (Intercept) 0.0956 0.3092
Number of obs: 100, groups: bloco, 5
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.09194 0.16634 0.553 0.5805
tratamentoF25 -0.22197 0.13080 -1.697 0.0897 .
tratamentoF50 -0.09389 0.13068 -0.719 0.4724
tratamentoN25 -0.54411 0.13224 -4.115 3.88e-05 ***
tratamentoN50 -0.17069 0.13072 -1.306 0.1916
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtF25 trtF50 trtN25
tratamntF25 -0.393
tratamntF50 -0.393 0.500
tratamntN25 -0.389 0.494 0.495
tratamntN50 -0.393 0.500 0.500 0.495
exp(-0.544)[1] 0.5804219
(52.3/(100-52.3))/ (39.2/(100-39.2))[1] 1.700595
Planted trees in treatment F25 have 19.91% less likely to survive (OR=exp(-0.222)) than in treatment 100% (\(p<0.1\))
Planted trees in treatment N25 have 41.96% less likely to survive (OR=exp(-0.544)) than in treatment 100% (\(p<0.001\))
round((1 - exp(-0.544))*100,2)[1] 41.96
emmeans::emmeans(mod_surv1, ~ tratamento, type = "response") |>
multcomp::cld(details = TRUE, adjust="sidak")$emmeans
tratamento prob SE df asymp.LCL asymp.UCL .group
N25 0.389 0.0398 Inf 0.293 0.495 1
F25 0.468 0.0414 Inf 0.364 0.574 12
N50 0.480 0.0415 Inf 0.376 0.586 2
F50 0.500 0.0416 Inf 0.394 0.605 2
100 0.523 0.0415 Inf 0.417 0.627 2
Confidence level used: 0.95
Conf-level adjustment: sidak method for 5 estimates
Intervals are back-transformed from the logit scale
P value adjustment: sidak method for 10 tests
Tests are performed on the log odds ratio scale
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
$comparisons
contrast odds.ratio SE df null z.ratio p.value
F25 / N25 1.38 0.183 Inf 1 2.436 0.1390
N50 / N25 1.45 0.192 Inf 1 2.825 0.0463
N50 / F25 1.05 0.138 Inf 1 0.392 1.0000
F50 / N25 1.57 0.207 Inf 1 3.407 0.0066
F50 / F25 1.14 0.149 Inf 1 0.980 0.9810
F50 / N50 1.08 0.141 Inf 1 0.588 0.9997
100 / N25 1.72 0.228 Inf 1 4.115 0.0004
100 / F25 1.25 0.163 Inf 1 1.697 0.6093
100 / N50 1.19 0.155 Inf 1 1.306 0.8809
100 / F50 1.10 0.144 Inf 1 0.719 0.9983
P value adjustment: sidak method for 10 tests
Tests are performed on the log odds ratio scale
survival_planted2024 |>
dplyr::group_by( tratamento)|>
dplyr::summarise(Mean = round(mean((n/24)*100),1),
Max = round(max((n/24)*100),1),
Min = round(min((n/24)*100),1) ) |>
dplyr::mutate(values = paste0(Mean," (",Min," - ",Max,")")) |>
dplyr::select(-c(Mean, Max, Min)) # A tibble: 5 × 2
tratamento values
<chr> <chr>
1 100 52.3 (29.2 - 79.2)
2 F25 46.9 (12.5 - 70.8)
3 F50 50 (29.2 - 75)
4 N25 39.2 (8.3 - 79.2)
5 N50 48.1 (20.8 - 75)
How planting scheme affected grass cover inside and outside planted areas?
Excluding treatment 100:
ordem <- c("N25", "F25","N50", "F50", "100")
grass_cover <- readxl::read_xlsx(
paste0(path, "/Itatinga_Nuc_grasscover.xlsx"))|>
dplyr::mutate(
grass_perc = dplyr::if_else(gramineas_perc>=1, 0.9999, gramineas_perc),
tratamento = ordered(tratamento,
levels = ordem,
labels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%"))) |>
droplevels()
grass_cover2024 <- grass_cover |> dplyr::filter(tratamento!="Planting 100%") |> dplyr::filter(ano==2024)
grass_cover2024 |>
with(table(bloco, parcela, by=tratamento)), , by = Nuclei 25%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Strips 25%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Nuclei 50%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Strips 50%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Planting 100%
parcela
bloco 1 2 3 4
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
Too many zeros? Zero-Inflated model? There is few zeros. Only binomial fit
grass_cover2024 |>
ggplot(aes(x=grass_perc))+
geom_histogram()+
facet_grid(plantio~tratamento)+theme_classic()`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
grass_cover2024 |>
with(table(bloco, parcela, by=tratamento)), , by = Nuclei 25%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Strips 25%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Nuclei 50%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Strips 50%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Planting 100%
parcela
bloco 1 2 3 4
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
grass_cover2024 <-grass_cover2024 |>
dplyr::mutate(Tratamento=factor(tratamento, ordered = FALSE))
modelo_gc24a <- glmmTMB::glmmTMB(grass_perc ~ Tratamento*plantio + (1 | bloco),
family = glmmTMB::beta_family(), data = grass_cover2024)
modelo_gc24b <- grass_cover2024 |>
dplyr::mutate(tratamento=factor(tratamento,ordered = FALSE))|>
with(glmmTMB::glmmTMB(grass_perc ~ tratamento+plantio + (1 | bloco),
ziformula = ~ 1, family = glmmTMB::beta_family()))Warning in glmmTMB::glmmTMB(grass_perc ~ tratamento + plantio + (1 | bloco), :
use of the 'data' argument is recommended
modelo_gc24c <- grass_cover2024 |>
dplyr::mutate(tratamento=factor(tratamento,ordered = FALSE))|>
with(glmmTMB::glmmTMB(grass_perc ~ plantio + (1 | bloco),
family = glmmTMB::beta_family()))Warning in glmmTMB::glmmTMB(grass_perc ~ plantio + (1 | bloco), family =
glmmTMB::beta_family()): use of the 'data' argument is recommended
modelo_gc24d <- grass_cover2024 |>
dplyr::mutate(tratamento=factor(tratamento,ordered = FALSE))|>
with(glmmTMB::glmmTMB(grass_perc ~ tratamento + (1 | bloco),
family = glmmTMB::beta_family()))Warning in glmmTMB::glmmTMB(grass_perc ~ tratamento + (1 | bloco), family =
glmmTMB::beta_family()): use of the 'data' argument is recommended
modelo_gc24null <- grass_cover2024 |>
dplyr::mutate(tratamento=factor(tratamento,ordered = FALSE))|>
with(glmmTMB::glmmTMB(grass_perc ~ 0 + (1 | bloco),
family = glmmTMB::beta_family()))Warning in glmmTMB::glmmTMB(grass_perc ~ 0 + (1 | bloco), family =
glmmTMB::beta_family()): use of the 'data' argument is recommended
AIC(modelo_gc24a, modelo_gc24b, modelo_gc24c, modelo_gc24d, modelo_gc24null) df AIC
modelo_gc24a 10 -61.25388
modelo_gc24b 8 -53.84400
modelo_gc24c 4 -51.58877
modelo_gc24d 6 -21.66566
modelo_gc24null 2 -10.55579
Analysis of residuals
par(mfrow=c(1,2))
DHARMa::plotQQunif(modelo_gc24a)
DHARMa::plotResiduals(modelo_gc24a)Is there any difference between treatments considering inside/outside?
resultado_24 <- emmeans::emmeans(modelo_gc24a, ~ Tratamento|plantio) |>
multcomp::cld( details = TRUE)
resultado_24$emmeansplantio = CP:
Tratamento emmean SE df asymp.LCL asymp.UCL .group
Strips 50% -0.727 0.228 Inf -1.173 -0.280 1
Strips 25% -0.300 0.224 Inf -0.739 0.138 12
Nuclei 25% 0.138 0.223 Inf -0.299 0.575 23
Nuclei 50% 0.708 0.228 Inf 0.262 1.154 3
plantio = SP:
Tratamento emmean SE df asymp.LCL asymp.UCL .group
Strips 50% 0.937 0.231 Inf 0.485 1.389 1
Nuclei 50% 0.945 0.231 Inf 0.493 1.397 1
Nuclei 25% 1.036 0.232 Inf 0.581 1.490 1
Strips 25% 1.166 0.234 Inf 0.708 1.625 1
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Results are given on the log odds ratio (not the response) scale.
P value adjustment: tukey method for comparing a family of 4 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Answer: yes.
Is there any difference between inside/outside considering treatments?
resultado_24a <- emmeans::emmeans(modelo_gc24a, ~ plantio|Tratamento) |>
multcomp::cld(details = TRUE)
resultado_24a$emmeansTratamento = Nuclei 25%:
plantio emmean SE df asymp.LCL asymp.UCL .group
CP 0.138 0.223 Inf -0.299 0.575 1
SP 1.036 0.232 Inf 0.581 1.490 2
Tratamento = Strips 25%:
plantio emmean SE df asymp.LCL asymp.UCL .group
CP -0.300 0.224 Inf -0.739 0.138 1
SP 1.166 0.234 Inf 0.708 1.625 2
Tratamento = Nuclei 50%:
plantio emmean SE df asymp.LCL asymp.UCL .group
CP 0.708 0.228 Inf 0.262 1.154 1
SP 0.945 0.231 Inf 0.493 1.397 1
Tratamento = Strips 50%:
plantio emmean SE df asymp.LCL asymp.UCL .group
CP -0.727 0.228 Inf -1.173 -0.280 1
SP 0.937 0.231 Inf 0.485 1.389 2
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Results are given on the log odds ratio (not the response) scale.
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Answer: yes.
broom.mixed::tidy(modelo_gc24a) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))|> knitr::kable()| effect | component | group | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 0.138 | 0.223 | 0.620 | 0.535 |
| fixed | cond | NA | TratamentoStrips 25% | -0.439 | 0.316 | -1.388 | 0.165 |
| fixed | cond | NA | TratamentoNuclei 50% | 0.569 | 0.318 | 1.790 | 0.073 |
| fixed | cond | NA | TratamentoStrips 50% | -0.865 | 0.319 | -2.710 | 0.007 |
| fixed | cond | NA | plantioSP | 0.897 | 0.321 | 2.794 | 0.005 |
| fixed | cond | NA | TratamentoStrips 25%:plantioSP | 0.569 | 0.453 | 1.256 | 0.209 |
| fixed | cond | NA | TratamentoNuclei 50%:plantioSP | -0.660 | 0.454 | -1.454 | 0.146 |
| fixed | cond | NA | TratamentoStrips 50%:plantioSP | 0.766 | 0.454 | 1.688 | 0.091 |
| ran_pars | cond | bloco | sd__(Intercept) | 0.000 | NA | NA | NA |
sjPlot::tab_model(modelo_gc24a)| grass perc | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 1.15 | 0.74 – 1.78 | 0.535 |
| Tratamento [Strips 25%] | 0.64 | 0.35 – 1.20 | 0.165 |
| Tratamento [Nuclei 50%] | 1.77 | 0.95 – 3.30 | 0.073 |
| Tratamento [Strips 50%] | 0.42 | 0.23 – 0.79 | 0.007 |
| plantio [SP] | 2.45 | 1.31 – 4.60 | 0.005 |
| Tratamento [Strips 25%] × plantio [SP] |
1.77 | 0.73 – 4.30 | 0.209 |
| Tratamento [Nuclei 50%] × plantio [SP] |
0.52 | 0.21 – 1.26 | 0.146 |
| Tratamento [Strips 50%] × plantio [SP] |
2.15 | 0.88 – 5.24 | 0.091 |
| Random Effects | |||
| σ2 | 0.16 | ||
| τ00 bloco | 0.00 | ||
| N bloco | 5 | ||
| Observations | 160 | ||
| Marginal R2 / Conditional R2 | 0.736 / NA | ||
lme4::VarCorr(modelo_gc24a) # ou summary(modelo_gc24a)$varcor
Conditional model:
Groups Name Std.Dev.
bloco (Intercept) 1.0464e-05
MuMIn::r.squaredGLMM(modelo_gc24a)Warning: Can't compute random effect variances. Some variance components equal
zero. Your model may suffer from singularity (see `?lme4::isSingular`
and `?performance::check_singularity`).
Decrease the `tolerance` level to force the calculation of random effect
variances, or impose priors on your random effects parameters (using
packages like `brms` or `glmmTMB`).
R2m R2c
[1,] 0.7355586 0.7355586
m_fix <- update(modelo_gc24a, . ~ . - (1|bloco))
anova(modelo_gc24a, m_fix) Data: grass_cover2024
Models:
m_fix: grass_perc ~ Tratamento + plantio + Tratamento:plantio, zi=~0, disp=~1
modelo_gc24a: grass_perc ~ Tratamento * plantio + (1 | bloco), zi=~0, disp=~1
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
m_fix 9 -63.254 -35.577 40.627 -81.254
modelo_gc24a 10 -61.254 -30.502 40.627 -81.254 0 1 1
performance::check_singularity(modelo_gc24a)[1] TRUE
Small letters: Comparison between treatments, only inside planted areas.
Capital letters: Comparison between treatments, only outside planted areas.
letras <- data.frame(
tratamento = rep(levels(as.factor(grass_cover$tratamento)),each=2),
plantio = levels(as.factor(grass_cover$plantio)),
label = c("BC ", " a",
"AB ", " a",
"C ", " a",
"A ", " a",
"A",""),
y = c(0.6, 0.82,
0.5, 0.8,
0.65, 0.8,
0.35, 0.8,
0.35, .1)
) |>
dplyr::mutate(
tratamento = forcats::fct_recode(tratamento,
"Nuclei\n25%" = "Nuclei 25%",
"Strips\n25%" = "Strips 25%",
"Nuclei\n50%" = "Nuclei 50%",
"Strips\n50%" = "Strips 50%",
"Planting\n100%" = "Planting 100%"
)
)
fig_grass <- grass_cover |>
dplyr::filter(ano==2024) |>
dplyr::group_by(tratamento, plantio) |>
dplyr::summarise(grass_cover_mean = mean(gramineas_perc),
sd = sd(gramineas_perc),
n = dplyr::n(),
se = sd/sqrt(n),
.groups="drop") |>
dplyr::mutate(
tratamento = forcats::fct_recode(tratamento,
"Nuclei\n25%" = "Nuclei 25%",
"Strips\n25%" = "Strips 25%",
"Nuclei\n50%" = "Nuclei 50%",
"Strips\n50%" = "Strips 50%",
"Planting\n100%" = "Planting 100%"
)
) |>
ggplot(aes(x = tratamento, y = grass_cover_mean, fill = plantio)) +
geom_bar(position = position_dodge2(preserve = "single"),
stat = "identity", color = "black") +
geom_errorbar(aes(ymin = grass_cover_mean - se, ymax = grass_cover_mean + se),
width = .2, position = position_dodge(0.9)) +
scale_fill_manual(
labels = c('Planted', 'Unplanted'),
values = c("#92D050", "white")) +
geom_text(data = letras, aes(x = tratamento, y = y, label = label), size=4, vjust = -0.5) +
labs(fill = "",
y = "Grass cover (%)",
x = "") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent_format(accuracy = 1),
limits = c(0,1)) +
theme_classic()+
theme(
axis.text.x = element_text(color="black", hjust = 0.5, size=11),
axis.text.y = element_text(color="black", size=11),
axis.title.y = element_text(color="black", size=13),
legend.position = c(0.84, 0.99),
legend.text = element_text(size = 11, margin = margin(0,0,0, l=1, "mm")),
legend.key.size = unit(3, "mm"),
legend.spacing = unit(0, "mm"),
legend.background = element_blank())Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.
fig_grass# ggsave("fig_grass.png", plot=fig_grass, width=90, height = 90 ,unit="mm")grass_cover |>
dplyr::group_by(ano, tratamento, plantio)|>
dplyr::summarise(Mean = mean(gramineas_perc)*100,
Max = max(gramineas_perc)*100,
Min = min(gramineas_perc)*100) |>
dplyr::mutate(values = paste0(Mean," (",Min," - ",Max,")")) |>
dplyr::select(-c(Mean, Max, Min)) |>
tidyr::pivot_wider(names_from = tratamento, values_from = values)`summarise()` has grouped output by 'ano', 'tratamento'. You can override using
the `.groups` argument.
# A tibble: 4 × 7
# Groups: ano [2]
ano plantio `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%`
<dbl> <chr> <chr> <chr> <chr> <chr>
1 2022 CP 29.2 (0 - 85) 16.9 (0 - 100) 34.9 (0 - 100) 12.35 (0 - …
2 2022 SP 65.25 (20 - 100) 65.175 (23 - 100) 51.6 (0 - 100) 60.35 (0 - …
3 2024 CP 54.07 (8.2 - 95) 39.91 (4 - 95) 56.3485 (2.64… 29.9795 (5 …
4 2024 SP 76.835 (39 - 97.3) 75.45 (29 - 100) 73.627 (26 - … 77.9725 (49…
# ℹ 1 more variable: `Planting 100%` <chr>
Only recruits (smaller and larger), not including planted individuals
species_table <- readxl::read_xlsx("itatinga_Nuc_regen_species_abundance.xlsx")
species_table|>
dplyr::select(Treatment, Block, Plot, Planted)|>
unique()|>
with(table(Block, Plot, by=Treatment)), , by = Full_Planted
Plot
Block P1 P2 P3 P4
B1 1 1 1 1
B2 1 1 1 1
B3 1 1 1 1
B4 1 1 1 1
B5 1 1 1 1
, , by = Nuclei_25
Plot
Block P1 P2 P3 P4
B1 2 2 2 2
B2 2 2 2 2
B3 2 2 2 2
B4 2 2 2 2
B5 2 2 2 2
, , by = Nuclei_50
Plot
Block P1 P2 P3 P4
B1 2 2 2 2
B2 2 2 2 2
B3 2 2 2 2
B4 2 2 2 2
B5 2 2 2 2
, , by = Strip_25
Plot
Block P1 P2 P3 P4
B1 2 2 2 2
B2 2 2 2 2
B3 2 2 2 2
B4 2 2 2 2
B5 2 2 2 2
, , by = Strip_50
Plot
Block P1 P2 P3 P4
B1 2 2 2 2
B2 2 2 2 2
B3 2 2 2 2
B4 2 2 2 2
B5 2 2 2 2
Viewing top species
total_smaller <- sum(species_table$n_sapling)
total_larger <- sum(species_table$n_tree)
species_table |>
dplyr::group_by(scientificName)|>
dplyr::summarise(
sum_large = sum(n_tree, na.rm=T)/total_larger,
sum_small = sum(n_sapling, na.rm=T)/total_smaller,
)|>
dplyr::arrange(desc(sum_small)) |>
head()# A tibble: 6 × 3
scientificName sum_large sum_small
<chr> <dbl> <dbl>
1 Baccharis dracunculifolia 0.0740 0.489
2 Vernonanthura polyanthes 0 0.313
3 Solanum sisymbriifolium 0.00157 0.0388
4 Croton floribundus 0.00157 0.0384
5 Senegalia polyphylla 0.00472 0.0305
6 Solanum granuloso-leprosum 0.0205 0.0233
Only small regenerants
dados_list <- species_table |>
dplyr::group_by(Treatment, Planted, scientificName) |>
dplyr::summarise(n_total = sum(n_sapling), .groups = "drop") |>
dplyr::group_by(Treatment, Planted) |>
dplyr::summarise(abundancias = list(n_total), .groups = "drop") |>
dplyr::mutate(nome_grupo = paste(Treatment, Planted, sep = "_")) %>%
{ setNames(.$abundancias, .$nome_grupo) }
res_inext <- iNEXT::iNEXT(
dados_list,
q = 0,
size=seq(0, 2000, 100),
datatype = "abundance"
)
estimativas <- res_inext$iNextEst |>
dplyr::bind_rows() |>
tidyr::separate(Assemblage, into = c("Treatment", "Percent", "Planted"), sep = "_")Figure Recruits’ species richness
For Hill q = 0
Only small regenerants and considering ruderal species (Vernonanthura & Baccharis)
dados_list <- species_table |>
dplyr::group_by(Treatment, Planted, scientificName) |>
dplyr::summarise(n_total = sum(n_sapling), .groups = "drop") |>
dplyr::group_by(Treatment, Planted) |>
dplyr::summarise(abundancias = list(n_total), .groups = "drop") |>
dplyr::mutate(nome_grupo = paste(Treatment, Planted, sep = "_")) %>%
{ setNames(.$abundancias, .$nome_grupo) }
res_inext <- iNEXT::iNEXT(
dados_list,
q = 0,
size=seq(0, 2000, 100),
datatype = "abundance"
)
estimativas <- res_inext$iNextEst |>
dplyr::bind_rows() |>
tidyr::separate(Assemblage, into = c("Treatment", "Percent", "Planted"), sep = "_")
curvas_df0 <- estimativas |>
dplyr::filter(Order.q==0)|>
dplyr::filter(!is.na(SC.LCL))|>
dplyr::mutate(
Percent = ifelse(Percent=="Planted", "100", Percent),
Percent = ordered(Percent, levels = c("25", "50", "100")),
Planted = ifelse(Planted=="Inside", "Planted", "Unplanted"),
Treatment = ifelse(Treatment=="Full", "Planting 100%", Treatment),
Treatment = factor(Treatment, levels = c("Nuclei", "Strip", "Planting 100%"), ordered = TRUE)
)
fig_recruit_rich <- curvas_df0 |>
ggplot(aes(x = m, y = qD, color = Treatment, linetype = Planted)) +
geom_line(linewidth = 1) +
geom_ribbon(aes(ymin = qD.LCL, ymax = qD.UCL, fill = Treatment), alpha = 0.2, colour = NA) +
geom_point(
data = dplyr::filter(curvas_df0, Method == "Observed"),
size = 2.5
) +
facet_grid(Percent~.,
labeller = as_labeller(
c("25" = "(a) Planting 25%",
"50" = "(b) Planting 50%",
"100" = "(c) Planting 100%")
)) +
scale_color_manual(values=c("red", "blue", "darkgreen"))+
scale_fill_manual(values=c("red", "blue", "darkgreen"))+
labs(
x = "Number of sampling units (regenerants)",
y = "Species diversity (Hill Index, q=0)"
) +
ggthemes::theme_clean()+
scale_x_continuous(expand = c(0.001, 0), breaks = seq(0, 3500, 500)) +
scale_y_continuous(expand = c(0, 0))+
theme(
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_blank(),
axis.text.x = element_text(color = "black", hjust = 1, size = 10),
axis.text.y = element_text(color = "black", size = 9),
axis.title.y = element_text(color = "black", size = 11),
strip.text = element_text(color = "black", hjust = 0.5, size = 10),
strip.placement = "outside",
legend.text = element_text(size = 11),
legend.title = element_blank(),
legend.key.size = unit(5, "mm"),
legend.box.background = element_blank(),
legend.box.margin = margin(t = 5, r = 0, b = 0, l = 0, "mm"),
legend.background = element_blank(),
legend.margin = margin(0, 0, 0, 0),
legend.position = c(0.25, 0.25),
panel.border = element_rect(colour = "black", fill = NA),
panel.spacing = unit(0, "lines")
)
fig_recruit_rich# ggsave("fig_richness_all_small_recruits.png", plot=fig_recruit_rich, width=90, height = 210 ,unit="mm")Estimating Density of recruits (\(ind.ha^{-1}\)) and considering ruderal species (Baccharis and Vernonanthura)
density_recruits_n_m2 <- species_table |>
dplyr::mutate(
n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
dplyr::group_by(Treatment, Block, Planted) |>
dplyr::summarise(
n_small_m2 = (sum(n_small, na.rm=T))/(144*4),
n_large_m2 = (sum(n_large, na.rm=T))/(144*4),
.groups = "drop")
dens_ha <- density_recruits_n_m2 |>
dplyr::mutate(
n_small_ha = n_small_m2 * 10000,
n_large_ha = n_large_m2 * 10000
)
dens_P_UP <- dens_ha |>
tidyr::pivot_wider(
names_from = Planted,
values_from = c(n_small_ha, n_large_ha),
names_sep = "_"
)
prop_table <- tidyr::tibble(
Treatment = unique(sort(species_table$Treatment)),
prop_P = c(1, 0.25, 0.5, 0.25, 0.5),
prop_UP = 1 - prop_P
)
result_regenerant <- dens_P_UP |>
dplyr::left_join(prop_table, by = "Treatment") |>
dplyr::mutate(
n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0) + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
) |>
dplyr::select(Treatment,
starts_with("n_small_ha"),
starts_with("n_large_ha"),
n_small_total_ha,
n_large_total_ha
)Smaller recruits
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_ha_Inside, na.rm=TRUE),1),
min = round(min(n_small_ha_Inside, na.rm=TRUE),0),
max = round(max(n_small_ha_Inside, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 3048.6 (2274.0 - 4705.0) 4871.5 (3490.0 - 7101.0) 6107.6 (2… 4816.0 … 6024.3 …
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
min = round(min(n_small_ha_Outside, na.rm=TRUE), 0),
max = round(max(n_small_ha_Outside, na.rm=TRUE), 0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_small_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 NaN (Inf - -Inf) 3961.8 (2344.0 - 8247.0) 3461.8 (2031.0 - … 3392.4 … 3871.5 …
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_total_ha, na.rm=TRUE),1),
min = round(min(n_small_total_ha, na.rm=TRUE),0),
max = round(max(n_small_total_ha, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 3048.6 (2274.0 - 4705.0) 2094.6 (872.0 - 6185.0) 2392.4 (10… 1874.1 … 2474.0 …
Larger recruits
All treatments and inside/outside areas showed similar values
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_ha_Inside, na.rm=TRUE),1),
min = round(min(n_large_ha_Inside, na.rm=TRUE),0),
max = round(max(n_large_ha_Inside, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 222.2 (104.0 - 486.0) 229.2 (0.0 - 608.0) 163.2 (35.0 - 295… 152.8 (… 225.7 (…
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_ha_Outside, na.rm=TRUE),1),
min = round(min(n_large_ha_Outside, na.rm=TRUE),0),
max = round(max(n_large_ha_Outside, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_large_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 NaN (Inf - -Inf) 364.6 (17.0 - 851.0) 163.2 (122.0 - 208.0) 302.1 (1… 381.9 (…
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_total_ha, na.rm = TRUE), 1),
min = round(min(n_large_total_ha, na.rm = TRUE), 0),
max = round(max(n_large_total_ha, na.rm = TRUE), 0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 222.2 (104.0 - 486.0) 165.4 (0.0 - 638.0) 81.6 (17.0 - 148.… 132.4 (… 151.9 (…
Data:
small_df <- species_table |>
dplyr::group_by(Treatment, Block, Plot, Planted)|>
dplyr::summarise(n = sum(n_sapling), .groups="drop")
small_df |>
ggplot(aes(x=Planted, y=n, color=Block))+
geom_point()+facet_wrap(~Treatment)+theme_bw()Models:
modelo_reg1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
# family = poisson(),
family = glmmTMB::nbinom2(),
data = small_df)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
modelo_reg2 <- glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
# family = poisson(),
family = glmmTMB::nbinom2(),
data = small_df)
modelo_reg3 <- glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
# family = poisson(),
family = glmmTMB::nbinom2(),
data = small_df)
modelo_reg_null <- glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
# family = poisson(), # Poisson
family = glmmTMB::nbinom2(), # NB
data = small_df)
# NBinom 1782.688 (null) | 1734.985 (Treatment*Planted) | 1756.665 (only Treament) | 1743.097 (only Planted)
# Poisson 4032.177 | 3373.854 | 3813.769 | 3759.763
AIC(modelo_reg_null, modelo_reg1, modelo_reg2, modelo_reg3) df AIC
modelo_reg_null 2 1782.688
modelo_reg1 11 1734.985
modelo_reg2 7 1756.665
modelo_reg3 4 1743.097
non-positive-definite Hessian matrix indicates that there are problems in estimating the model parameters. This may be related to redundancies in the random effects, perfect separation, few data for some combinations of factors, or problems in choosing the distribution. Therefore, let’s check the singularity:
performance::check_singularity(modelo_reg1)[1] FALSE
performance::check_model(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
`check_outliers()` does not yet support models of class `glmmTMB`.
Exploring residuals
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
DHARMa::plotQQunif(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
DHARMa::plotResiduals(modelo_reg3)
DHARMa::plotQQunif(modelo_reg3)car::Anova(modelo_reg1)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: n
Chisq Df Pr(>Chisq)
Treatment 23.4883 4 0.0001011 ***
Planted 30.8344 1 2.81e-08 ***
Treatment:Planted 1.8521 3 0.6036618
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Is there any difference between treatments considering only inside or outside?
emmeans::emmeans(modelo_reg1, ~ Treatment|Planted) |>
multcomp::cld(adjust="bonferroni")Planted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Full_Planted 3.76 0.175 Inf 3.31 4.21 1
Strip_25 4.24 0.173 Inf 3.79 4.68 2
Nuclei_25 4.24 0.174 Inf 3.80 4.69 2
Nuclei_50 4.40 0.173 Inf 3.96 4.85 2
Strip_50 4.42 0.173 Inf 3.98 4.87 2
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Strip_25 3.82 0.175 Inf 3.38 4.25 1
Nuclei_50 3.82 0.175 Inf 3.39 4.26 1
Nuclei_25 3.96 0.174 Inf 3.52 4.39 1
Strip_50 3.98 0.174 Inf 3.55 4.42 1
Full_Planted nonEst NA NA NA NA
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
P value adjustment: bonferroni method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Answer: Yes (inside) and no (outside).
Is there any difference between treatments not considering inside/outside?
emmeans::emmeans(modelo_reg2, ~ Treatment) |>
multcomp::cld(adjust="bonferroni") Treatment emmean SE df asymp.LCL asymp.UCL .group
Full_Planted 3.76 0.177 Inf 3.30 4.22 1
Strip_25 4.05 0.155 Inf 3.65 4.45 12
Nuclei_25 4.11 0.155 Inf 3.71 4.51 12
Nuclei_50 4.16 0.155 Inf 3.76 4.56 12
Strip_50 4.23 0.154 Inf 3.83 4.63 2
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 5 estimates
P value adjustment: bonferroni method for 10 tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Answer: No.
Is there any difference between inside/outside considering treatments?
emmeans::emmeans(modelo_reg1, ~ Planted|Treatment) |>
multcomp::cld(adjust="bonferroni")Treatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 3.76 0.175 Inf 3.42 4.10 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.96 0.174 Inf 3.57 4.35 1
Inside 4.24 0.174 Inf 3.85 4.63 1
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.82 0.175 Inf 3.43 4.21 1
Inside 4.40 0.173 Inf 4.02 4.79 2
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.82 0.175 Inf 3.42 4.21 1
Inside 4.24 0.173 Inf 3.85 4.63 2
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.98 0.174 Inf 3.59 4.37 1
Inside 4.42 0.173 Inf 4.03 4.81 2
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Answer: Yes.
sjPlot::tab_model(modelo_reg1)Model matrix is rank deficient. Parameters
`TreatmentStrip_50:PlantedOutside` were not estimable.
| n | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 42.99 | 30.52 – 60.55 | <0.001 |
| Treatment [Nuclei_25] | 1.62 | 1.19 – 2.20 | 0.002 |
| Treatment [Nuclei_50] | 1.90 | 1.40 – 2.59 | <0.001 |
| Treatment [Strip_25] | 1.61 | 1.19 – 2.19 | 0.002 |
| Treatment [Strip_50] | 1.94 | 1.43 – 2.63 | <0.001 |
| Planted [Outside] | 0.64 | 0.48 – 0.87 | 0.005 |
| Treatment [Nuclei_25] × Planted [Outside] |
1.17 | 0.76 – 1.80 | 0.478 |
| Treatment [Nuclei_50] × Planted [Outside] |
0.87 | 0.56 – 1.33 | 0.515 |
| Treatment [Strip_25] × Planted [Outside] |
1.02 | 0.66 – 1.57 | 0.937 |
| Random Effects | |||
| σ2 | 0.21 | ||
| τ00 Block | 0.09 | ||
| ICC | 0.30 | ||
| N Block | 5 | ||
| Observations | 180 | ||
| Marginal R2 / Conditional R2 | 0.164 / 0.412 | ||
broom.mixed::tidy(modelo_reg1) |>
dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |> dplyr::select(-c(group,component))# A tibble: 11 × 6
effect term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 fixed (Intercept) 3.76 0.175 21.5 0
2 fixed TreatmentNuclei_25 0.481 0.156 3.08 0.002
3 fixed TreatmentNuclei_50 0.644 0.156 4.13 0
4 fixed TreatmentStrip_25 0.476 0.156 3.05 0.002
5 fixed TreatmentStrip_50 0.662 0.156 4.25 0
6 fixed PlantedOutside -0.439 0.155 -2.83 0.005
7 fixed TreatmentNuclei_25:PlantedOuts… 0.156 0.219 0.709 0.478
8 fixed TreatmentNuclei_50:PlantedOuts… -0.143 0.22 -0.652 0.515
9 fixed TreatmentStrip_25:PlantedOutsi… 0.017 0.22 0.079 0.937
10 fixed TreatmentStrip_50:PlantedOutsi… NA NA NA NA
11 ran_pars sd__(Intercept) 0.301 NA NA NA
Data:
large_regenerants_2024 <- species_table |>
dplyr::group_by(Treatment, Block, Plot, Planted)|>
dplyr::summarise(n = sum(n_tree), .groups="drop")Plotting
large_regenerants_2024 |>
ggplot(aes(x=Planted, y=n, color=Block))+
geom_point()+facet_wrap(~Treatment)+theme_bw()Models
modelo_reg1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = large_regenerants_2024)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
modelo_reg2 <- glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = large_regenerants_2024)
modelo_reg3 <- glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = large_regenerants_2024)
modelo_reg_null <- glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
ziformula = ~1,
# family = poisson(), # ZIP
family = glmmTMB::nbinom2(), # ZINB
# family = glmmTMB::truncated_nbinom2(), # hurdle
data = large_regenerants_2024)
# Model (null) | (treat*planted) | (treat) | (planted)
# ZINB 885.72 | 869.63 | 868.93 | 864.21
# ZIP 1058.14 | 1015.06 | 1025.38 | 1024.49
# Hurdle 880.75 | 872.21 | 868.74 | 975.88
AIC(modelo_reg_null, modelo_reg1, modelo_reg2, modelo_reg3) df AIC
modelo_reg_null 3 885.7221
modelo_reg1 12 869.6264
modelo_reg2 8 868.9342
modelo_reg3 5 864.2072
performance::check_model(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
`check_outliers()` does not yet support models of class `glmmTMB`.
Exploring residuals
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
DHARMa::plotQQunif(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
DHARMa::plotResiduals(modelo_reg3)
DHARMa::plotQQunif(modelo_reg3)MuMIn::r.squaredGLMM(modelo_reg3)Warning: the null model is only correct if all the variables it uses are identical
to those used in fitting the original model.
R2m R2c
delta 0.03380077 0.03380078
lognormal 0.05135140 0.05135140
trigamma 0.01854201 0.01854201
# 0.03253385/0.036655970.0328/ 0.0338 variance explained by fixed effects and by the entire model (around 0.9999997 % is explained by fixed effects)
performance::r2(modelo_reg1)Warning: Can't compute random effect variances. Some variance components equal
zero. Your model may suffer from singularity (see `?lme4::isSingular`
and `?performance::check_singularity`).
Decrease the `tolerance` level to force the calculation of random effect
variances, or impose priors on your random effects parameters (using
packages like `brms` or `glmmTMB`).
[1] NA
broom.mixed::tidy(modelo_reg1) |>
dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |>
dplyr::select(-c(component, group))|>
knitr::kable()| effect | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|
| fixed | (Intercept) | 1.175 | 0.251 | 4.690 | 0.000 |
| fixed | TreatmentNuclei_25 | 0.064 | 0.369 | 0.173 | 0.863 |
| fixed | TreatmentNuclei_50 | -0.286 | 0.369 | -0.776 | 0.438 |
| fixed | TreatmentStrip_25 | -0.353 | 0.369 | -0.957 | 0.339 |
| fixed | TreatmentStrip_50 | 0.028 | 0.354 | 0.080 | 0.936 |
| fixed | PlantedOutside | 0.510 | 0.346 | 1.476 | 0.140 |
| fixed | TreatmentNuclei_25:PlantedOutside | -0.045 | 0.492 | -0.092 | 0.927 |
| fixed | TreatmentNuclei_50:PlantedOutside | -0.515 | 0.507 | -1.018 | 0.309 |
| fixed | TreatmentStrip_25:PlantedOutside | 0.169 | 0.499 | 0.339 | 0.734 |
| fixed | TreatmentStrip_50:PlantedOutside | NA | NA | NA | NA |
| fixed | (Intercept) | -3.510 | 3.051 | -1.150 | 0.250 |
| ran_pars | sd__(Intercept) | 0.000 | NA | NA | NA |
car::Anova(modelo_reg1)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: n
Chisq Df Pr(>Chisq)
Treatment 7.3216 4 0.11984
Planted 5.5546 1 0.01843 *
Treatment:Planted 1.9111 3 0.59106
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Is there any difference between treatments considering inside/outside?
emmeans::emmeans(modelo_reg1, ~ Treatment|Planted) |>
multcomp::cld(adjust = "bonferroni")Planted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Strip_25 0.822 0.284 Inf 0.0899 1.55 1
Nuclei_50 0.889 0.284 Inf 0.1589 1.62 1
Full_Planted 1.175 0.251 Inf 0.5298 1.82 1
Strip_50 1.203 0.260 Inf 0.5327 1.87 1
Nuclei_25 1.239 0.287 Inf 0.4988 1.98 1
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Nuclei_50 0.884 0.277 Inf 0.1920 1.58 1
Strip_25 1.501 0.259 Inf 0.8532 2.15 1
Nuclei_25 1.704 0.276 Inf 1.0145 2.39 1
Strip_50 1.714 0.235 Inf 1.1262 2.30 1
Full_Planted nonEst NA NA NA NA
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
P value adjustment: bonferroni method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Answer: No.
Is there any difference between inside/outside considering treatments?
emmeans::emmeans(modelo_reg1, ~ Planted|Treatment) |>
multcomp::cld(adjust = "bonferroni")Treatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.175 0.251 Inf 0.684 1.67 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.239 0.287 Inf 0.595 1.88 1
Outside 1.704 0.276 Inf 1.085 2.32 1
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 0.884 0.277 Inf 0.263 1.51 1
Inside 0.889 0.284 Inf 0.254 1.52 1
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 0.822 0.284 Inf 0.185 1.46 1
Outside 1.501 0.259 Inf 0.920 2.08 1
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.203 0.260 Inf 0.620 1.79 1
Outside 1.714 0.235 Inf 1.186 2.24 1
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Answer: No.
Getting proportions of planted/unplanted species
prop_regenerant_full <- species_table |>
dplyr::mutate(
planted_tree = ifelse(`Planting group`=="NA", "No", "Yes"))|>
dplyr::filter(n_sapling>0)|>
dplyr::group_by(Treatment, Planted, Block, Plot) |>
dplyr::summarise(
n_species = dplyr::n(),
n_ind = sum(n_sapling),
n_planted_ind = sum(ifelse(planted_tree == "Yes", n_sapling, 0), na.rm = TRUE),
n_planted_ind = ifelse(is.na(n_planted_ind), 0, n_planted_ind),
n_planted_spp = sum(planted_tree == "Yes"),
n_planted_spp = ifelse(is.na(n_planted_spp), 0, n_planted_spp),
.groups = "drop") |>
dplyr::mutate(
prop_planted_ind = n_planted_ind / n_ind,
prop_planted_spp = n_planted_spp / n_species
)
prop_regenerant_full |>
dplyr::group_by(Treatment, Planted)|>
dplyr::summarise(
prop_planted_mean = mean(prop_planted_ind, na.rm=T),
prop_planted_min = min(prop_planted_ind, na.rm=T),
prop_planted_max = max(prop_planted_ind, na.rm=T),
prop_planted_spp_mean = mean(prop_planted_spp, na.rm=T),
prop_planted_spp_min = min(prop_planted_spp, na.rm=T),
prop_planted_spp_max = max(prop_planted_spp, na.rm=T),
.groups = "drop"
) |>
dplyr::mutate(
prop_planted = sprintf("%.1f (%.1f - %.1f)",
round(prop_planted_mean,3)*100,
round(prop_planted_min,3)*100,
round(prop_planted_max,3)*100),
prop_planted_spp = sprintf("%.1f (%.1f - %.1f)",
round(prop_planted_spp_mean,3)*100,
round(prop_planted_spp_min,3)*100,
round(prop_planted_spp_max,3)*100)
) |>
dplyr::select(Treatment,Planted, prop_planted, prop_planted_spp) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = c(prop_planted, prop_planted_spp)
) # |> writexl::write_xlsx("resulta_prop_recruit.xlsx")# A tibble: 2 × 11
Planted prop_planted_Full_Plan…¹ prop_planted_Nuclei_25 prop_planted_Nuclei_50
<chr> <chr> <chr> <chr>
1 Inside 10.5 (0.0 - 43.3) 3.9 (0.0 - 14.8) 7.7 (0.0 - 44.4)
2 Outside <NA> 10.4 (0.0 - 91.2) 10.0 (0.0 - 69.2)
# ℹ abbreviated name: ¹prop_planted_Full_Planted
# ℹ 7 more variables: prop_planted_Strip_25 <chr>, prop_planted_Strip_50 <chr>,
# prop_planted_spp_Full_Planted <chr>, prop_planted_spp_Nuclei_25 <chr>,
# prop_planted_spp_Nuclei_50 <chr>, prop_planted_spp_Strip_25 <chr>,
# prop_planted_spp_Strip_50 <chr>
Statistical analysis
mod_binomial_full <-lme4::glmer(
cbind(n_planted_ind, n_ind - n_planted_ind) ~ Treatment*Planted + (1 | Block),
family = binomial(link = "logit"),
data = prop_regenerant_full)fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
performance::check_overdispersion(mod_binomial_full) # no# Overdispersion test
dispersion ratio = 1.748
p-value = 0.272
No overdispersion detected.
performance::check_model(mod_binomial_full)DHARMa::plotResiduals(mod_binomial_full)emmeans::emmeans(mod_binomial_full, pairwise ~ Treatment | Planted) |> multcomp::cld(adjust="sidak")Planted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Nuclei_25 -3.35 0.292 Inf -4.10 -2.598 1
Strip_50 -2.63 0.271 Inf -3.32 -1.933 2
Nuclei_50 -2.59 0.269 Inf -3.28 -1.901 2
Strip_25 -2.40 0.269 Inf -3.09 -1.708 2
Full_Planted -2.30 0.279 Inf -3.02 -1.583 2
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Strip_25 -2.79 0.285 Inf -3.50 -2.084 1
Nuclei_25 -2.48 0.278 Inf -3.17 -1.789 12
Nuclei_50 -2.30 0.276 Inf -2.99 -1.610 2
Strip_50 -1.56 0.265 Inf -2.22 -0.903 3
Full_Planted nonEst NA NA NA NA
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: sidak method for varying numbers of estimates
Results are given on the log odds ratio (not the response) scale.
P value adjustment: sidak method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
emmeans::emmeans(mod_binomial_full, pairwise ~ Planted | Treatment) |> multcomp::cld(adjust="sidak")Treatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside -2.30 0.279 Inf -2.85 -1.753 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside -3.35 0.292 Inf -4.00 -2.695 1
Outside -2.48 0.278 Inf -3.10 -1.860 2
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside -2.59 0.269 Inf -3.19 -1.991 1
Outside -2.30 0.276 Inf -2.92 -1.681 2
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -2.79 0.285 Inf -3.43 -2.156 1
Inside -2.40 0.269 Inf -3.00 -1.797 2
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside -2.63 0.271 Inf -3.23 -2.023 1
Outside -1.56 0.265 Inf -2.16 -0.971 2
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: sidak method for varying numbers of estimates
Results are given on the log odds ratio (not the response) scale.
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
broom.mixed::tidy(mod_binomial_full) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 10 × 7
effect group term estimate std.error statistic p.value
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) -2.30 0.279 -8.24 0
2 fixed <NA> TreatmentNuclei_25 -1.05 0.188 -5.56 0
3 fixed <NA> TreatmentNuclei_50 -0.293 0.151 -1.95 0.051
4 fixed <NA> TreatmentStrip_25 -0.101 0.151 -0.665 0.506
5 fixed <NA> TreatmentStrip_50 -0.329 0.154 -2.14 0.033
6 fixed <NA> PlantedOutside 1.06 0.128 8.33 0
7 fixed <NA> TreatmentNuclei_25:Plant… -0.199 0.226 -0.88 0.379
8 fixed <NA> TreatmentNuclei_50:Plant… -0.77 0.193 -3.99 0
9 fixed <NA> TreatmentStrip_25:Plante… -1.46 0.205 -7.12 0
10 ran_pars Block sd__(Intercept) 0.564 NA NA NA
sjPlot::tab_model(mod_binomial_full)| cbind(n planted ind,n ind-n planted ind) |
|||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.10 | 0.06 – 0.17 | <0.001 |
| Treatment [Nuclei_25] | 0.35 | 0.24 – 0.51 | <0.001 |
| Treatment [Nuclei_50] | 0.75 | 0.56 – 1.00 | 0.051 |
| Treatment [Strip_25] | 0.90 | 0.67 – 1.22 | 0.506 |
| Treatment [Strip_50] | 0.72 | 0.53 – 0.97 | 0.033 |
| Planted [Outside] | 2.90 | 2.26 – 3.72 | <0.001 |
| Treatment [Nuclei_25] × Planted [Outside] |
0.82 | 0.53 – 1.28 | 0.379 |
| Treatment [Nuclei_50] × Planted [Outside] |
0.46 | 0.32 – 0.68 | <0.001 |
| Treatment [Strip_25] × Planted [Outside] |
0.23 | 0.16 – 0.35 | <0.001 |
| Random Effects | |||
| σ2 | 0.05 | ||
| τ00 Block | 0.32 | ||
| ICC | 0.86 | ||
| N Block | 5 | ||
| Observations | 178 | ||
| Marginal R2 / Conditional R2 | 0.354 / 0.910 | ||
Hill q = 0
Only small regenerants and considering ruderal species (Vernonanthura & Baccharis)
dados_list <- species_table |>
dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
dplyr::group_by(Treatment, Planted, scientificName) |>
dplyr::summarise(n_total = sum(n_sapling), .groups = "drop") |>
dplyr::group_by(Treatment, Planted) |>
dplyr::summarise(abundancias = list(n_total), .groups = "drop") |>
dplyr::mutate(nome_grupo = paste(Treatment, Planted, sep = "_")) %>%
{ setNames(.$abundancias, .$nome_grupo) }
res_inext <- iNEXT::iNEXT(
dados_list,
q = 0,
size=seq(0, 450, 10),
datatype = "abundance"
)
estimativas <- res_inext$iNextEst |>
dplyr::bind_rows() |>
tidyr::separate(Assemblage, into = c("Treatment", "Percent", "Planted"), sep = "_")
curvas_df0 <- estimativas |>
dplyr::filter(Order.q==0)|>
dplyr::filter(!is.na(SC.LCL))|>
dplyr::mutate(
Percent = ifelse(Percent=="Planted", "100", Percent),
Percent = ordered(Percent, levels = c("25", "50", "100")),
Planted = ifelse(Planted=="Inside", "Planted", "Unplanted"),
Treatment = ifelse(Treatment=="Full", "Planting 100%", Treatment),
Treatment = factor(Treatment, levels = c("Nuclei", "Strip", "Planting 100%"), ordered = TRUE)
)
fig_recruit_rich <- curvas_df0 |>
ggplot(aes(x = m, y = qD, color = Treatment, linetype = Planted)) +
geom_line(linewidth = 1) +
geom_ribbon(aes(ymin = qD.LCL, ymax = qD.UCL, fill = Treatment), alpha = 0.2, colour = NA) +
geom_point(
data = dplyr::filter(curvas_df0, Method == "Observed"),
size = 2.5
) +
facet_grid(Percent~.,
labeller = as_labeller(
c("25" = "(a) Planting 25%",
"50" = "(b) Planting 50%",
"100" = "(c) Planting 100%")
)) +
scale_color_manual(values=c("red", "blue", "darkgreen"))+
scale_fill_manual(values=c("red", "blue", "darkgreen"))+
labs(
x = "Number of sampling units (regenerants)",
y = "Species diversity (Hill Index, q=0)"
) +
ggthemes::theme_clean()+
scale_x_continuous(expand = c(0.001, 0), breaks = seq(0, 450, 50)) +
scale_y_continuous(expand = c(0, 0))+
theme(
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_blank(),
axis.text.x = element_text(color = "black", hjust = 1, size = 10),
axis.text.y = element_text(color = "black", size = 9),
axis.title.y = element_text(color = "black", size = 11),
strip.text = element_text(color = "black", hjust = 0.5, size = 10),
strip.placement = "outside",
legend.text = element_text(size = 11),
legend.title = element_blank(),
legend.key.size = unit(5, "mm"),
legend.box.background = element_blank(),
legend.box.margin = margin(t = 5, r = 0, b = 0, l = 0, "mm"),
legend.background = element_blank(),
legend.margin = margin(0, 0, 0, 0),
legend.position = c(0.25, 0.25),
panel.border = element_rect(colour = "black", fill = NA),
panel.spacing = unit(0, "lines")
)Plotting and saving
fig_recruit_richggsave("fig_richness_no_ruderal_small_recruits.png", plot=fig_recruit_rich, width=90, height = 210 ,unit="mm")Estimating Density of recruits (\(ind.ha^{-1}\)) without Baccharis and Vernonanthura
density_recruits_n_m2 <- species_table |>
dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
dplyr::mutate(
n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
dplyr::group_by(Treatment, Block, Planted) |>
dplyr::summarise(
n_small_m2 = (sum(n_small, na.rm=T))/(144*4),
n_large_m2 = (sum(n_large, na.rm=T))/(144*4),
.groups = "drop")
dens_ha <- density_recruits_n_m2 |>
dplyr::mutate(
n_small_ha = n_small_m2 * 10000,
n_large_ha = n_large_m2 * 10000
)
dens_P_UP <- dens_ha |>
tidyr::pivot_wider(
names_from = Planted,
values_from = c(n_small_ha, n_large_ha),
names_sep = "_"
)
prop_table <- tidyr::tibble(
Treatment = unique(sort(species_table$Treatment)),
prop_P = c(1, 0.25, 0.5, 0.25, 0.5),
prop_UP = 1 - prop_P
)
result_regenerant <- dens_P_UP |>
dplyr::left_join(prop_table, by = "Treatment") |>
dplyr::mutate(
n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0) + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
) |>
dplyr::select(Treatment,
starts_with("n_small_ha"),
starts_with("n_large_ha"),
n_small_total_ha,
n_large_total_ha
)Smaller recruits
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_ha_Inside, na.rm=TRUE),1),
min = round(min(n_small_ha_Inside, na.rm=TRUE),0),
max = round(max(n_small_ha_Inside, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 503.5 (139.0 - 1146.0) 482.6 (174.0 - 903.0) 961.8 (278.0 -… 722.2 (… 614.6 (…
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
min = round(min(n_small_ha_Outside, na.rm=TRUE), 0),
max = round(max(n_small_ha_Outside, na.rm=TRUE), 0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_small_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 NaN (Inf - -Inf) 1177.1 (469.0 - 1910.0) 822.9 (538.0 - 138… 1145.8 … 1381.9 …
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_total_ha, na.rm=TRUE),1),
min = round(min(n_small_total_ha, na.rm=TRUE),0),
max = round(max(n_small_total_ha, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 503.5 (139.0 - 1146.0) 501.7 (43.0 - 1432.0) 446.2 (139.0 -… 577.7 (… 499.1 (…
Larger recruits
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_ha_Inside, na.rm=TRUE),1),
min = round(min(n_large_ha_Inside, na.rm=TRUE),0),
max = round(max(n_large_ha_Inside, na.rm=TRUE),0)
) |>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 208.3 (87.0 - 486.0) 211.8 (0.0 - 608.0) 159.7 (35.0 - 278.… 149.3 (… 222.2 (…
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_ha_Outside, na.rm=TRUE),1),
min = round(min(n_large_ha_Outside, na.rm=TRUE),0),
max = round(max(n_large_ha_Outside, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_large_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 NaN (Inf - -Inf) 333.3 (0.0 - 851.0) 138.9 (104.0 - 208.0) 253.5 (35… 364.6 (…
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_total_ha, na.rm = TRUE), 1),
min = round(min(n_large_total_ha, na.rm = TRUE), 0),
max = round(max(n_large_total_ha, na.rm = TRUE), 0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 208.3 (87.0 - 486.0) 151.5 (0.0 - 638.0) 74.7 (17.0 - 139.0) 126.4 (… 146.7 (…
Data:
small_df <- species_table |>
dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
dplyr::group_by(Treatment, Block, Plot, Planted)|>
dplyr::summarise(n = sum(n_sapling), .groups="drop")
small_df |>
ggplot(aes(x=Planted, y=n, color=Block))+
geom_point()+facet_wrap(~Treatment)+theme_bw()Models:
modelo_reg1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
# family = poisson(),
family = glmmTMB::nbinom2(),
data = small_df)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
modelo_reg2 <- glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
# family = poisson(),
family = glmmTMB::nbinom2(),
data = small_df)
modelo_reg3 <- glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
# family = poisson(),
family = glmmTMB::nbinom2(),
data = small_df)
modelo_reg_null <- glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
# family = poisson(), # Poisson
family = glmmTMB::nbinom2(), # NB
data = small_df)
# NBinom 1270.77 (null) | 1226.20 (Treatment*Planted) | 1250.26 (only Treament) | 1227.00 (only Planted)
# Poisson 4032.177 | 3373.854 | 3813.769 | 3759.763
AIC(modelo_reg_null, modelo_reg1, modelo_reg2, modelo_reg3) df AIC
modelo_reg_null 2 1270.773
modelo_reg1 11 1226.199
modelo_reg2 7 1250.263
modelo_reg3 4 1227.001
non-positive-definite Hessian matrix indicates that there are problems in estimating the model parameters. This may be related to redundancies in the random effects, perfect separation, few data for some combinations of factors, or problems in choosing the distribution. Therefore, let’s check the singularity:
performance::check_singularity(modelo_reg1)[1] FALSE
performance::check_model(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
`check_outliers()` does not yet support models of class `glmmTMB`.
Exploring residuals
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
DHARMa::plotQQunif(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
DHARMa::plotResiduals(modelo_reg3)
DHARMa::plotQQunif(modelo_reg3)car::Anova(modelo_reg1)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: n
Chisq Df Pr(>Chisq)
Treatment 4.7099 4 0.31838
Planted 24.8150 1 6.311e-07 ***
Treatment:Planted 10.7731 3 0.01302 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Is there any difference between treatments considering only inside or outside?
emmeans::emmeans(modelo_reg1, ~ Treatment|Planted) |>
multcomp::cld(adjust="bonferroni")Planted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Full_Planted 1.86 0.241 Inf 1.24 2.48 1
Nuclei_25 1.87 0.241 Inf 1.25 2.50 1
Strip_50 2.09 0.238 Inf 1.48 2.70 1
Strip_25 2.27 0.236 Inf 1.66 2.88 1
Nuclei_50 2.42 0.234 Inf 1.82 3.02 1
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Nuclei_50 2.43 0.234 Inf 1.84 3.01 1
Strip_25 2.73 0.232 Inf 2.15 3.31 1
Nuclei_25 2.79 0.231 Inf 2.21 3.36 1
Strip_50 2.99 0.230 Inf 2.41 3.56 1
Full_Planted nonEst NA NA NA NA
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
P value adjustment: bonferroni method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Answer: No
Is there any difference between inside/outside considering treatments?
emmeans::emmeans(modelo_reg1, ~ Planted|Treatment) |>
multcomp::cld(adjust="bonferroni")Treatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.86 0.241 Inf 1.39 2.33 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.87 0.241 Inf 1.33 2.42 1
Outside 2.79 0.231 Inf 2.27 3.30 2
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 2.42 0.234 Inf 1.90 2.95 1
Outside 2.43 0.234 Inf 1.90 2.95 1
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 2.27 0.236 Inf 1.74 2.80 1
Outside 2.73 0.232 Inf 2.21 3.25 2
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 2.09 0.238 Inf 1.56 2.63 1
Outside 2.99 0.230 Inf 2.47 3.51 2
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Answer: Yes.
sjPlot::tab_model(modelo_reg1)Model matrix is rank deficient. Parameters
`TreatmentStrip_50:PlantedOutside` were not estimable.
| n | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 6.43 | 4.01 – 10.31 | <0.001 |
| Treatment [Nuclei_25] | 1.01 | 0.63 – 1.62 | 0.955 |
| Treatment [Nuclei_50] | 1.75 | 1.11 – 2.76 | 0.016 |
| Treatment [Strip_25] | 1.51 | 0.95 – 2.39 | 0.079 |
| Treatment [Strip_50] | 1.26 | 0.79 – 2.00 | 0.326 |
| Planted [Outside] | 2.45 | 1.57 – 3.82 | <0.001 |
| Treatment [Nuclei_25] × Planted [Outside] |
1.01 | 0.54 – 1.91 | 0.964 |
| Treatment [Nuclei_50] × Planted [Outside] |
0.41 | 0.22 – 0.77 | 0.005 |
| Treatment [Strip_25] × Planted [Outside] |
0.65 | 0.35 – 1.21 | 0.171 |
| Random Effects | |||
| σ2 | 0.40 | ||
| τ00 Block | 0.15 | ||
| ICC | 0.27 | ||
| N Block | 5 | ||
| Observations | 180 | ||
| Marginal R2 / Conditional R2 | 0.207 / 0.419 | ||
broom.mixed::tidy(modelo_reg1) |>
dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |> dplyr::select(-c(group,component))# A tibble: 11 × 6
effect term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 fixed (Intercept) 1.86 0.241 7.72 0
2 fixed TreatmentNuclei_25 0.014 0.24 0.056 0.955
3 fixed TreatmentNuclei_50 0.56 0.232 2.41 0.016
4 fixed TreatmentStrip_25 0.411 0.234 1.75 0.079
5 fixed TreatmentStrip_50 0.232 0.236 0.982 0.326
6 fixed PlantedOutside 0.897 0.226 3.97 0
7 fixed TreatmentNuclei_25:PlantedOuts… 0.015 0.322 0.045 0.964
8 fixed TreatmentNuclei_50:PlantedOuts… -0.892 0.319 -2.80 0.005
9 fixed TreatmentStrip_25:PlantedOutsi… -0.437 0.319 -1.37 0.171
10 fixed TreatmentStrip_50:PlantedOutsi… NA NA NA NA
11 ran_pars sd__(Intercept) 0.383 NA NA NA
Data:
large_regenerants_2024 <- species_table |>
dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
dplyr::group_by(Treatment, Block, Plot, Planted)|>
dplyr::summarise(n = sum(n_tree), .groups="drop")Plotting
large_regenerants_2024 |>
ggplot(aes(x=Planted, y=n, color=Block))+
geom_point()+facet_wrap(~Treatment)+theme_bw()Models
modelo_reg1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = large_regenerants_2024)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
modelo_reg2 <- glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = large_regenerants_2024)
modelo_reg3 <- glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = large_regenerants_2024)
modelo_reg_null <- glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
ziformula = ~1,
# family = poisson(), # ZIP
family = glmmTMB::nbinom2(), # ZINB
# family = glmmTMB::truncated_nbinom2(), # hurdle
data = large_regenerants_2024)
# Model (null) | (treat*planted) | (treat) | (planted)
# ZINB 855.21 | 843.24 | 840.70 | 838.42
AIC(modelo_reg_null, modelo_reg1, modelo_reg2, modelo_reg3) df AIC
modelo_reg_null 3 855.2141
modelo_reg1 12 843.2460
modelo_reg2 8 840.7029
modelo_reg3 5 838.4200
performance::check_model(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
`check_outliers()` does not yet support models of class `glmmTMB`.
Exploring residuals
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
DHARMa::plotQQunif(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
DHARMa::plotResiduals(modelo_reg3)
DHARMa::plotQQunif(modelo_reg3)car::Anova(modelo_reg1)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: n
Chisq Df Pr(>Chisq)
Treatment 7.9554 4 0.09323 .
Planted 3.4787 1 0.06216 .
Treatment:Planted 2.1403 3 0.54381
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Is there any difference between treatments considering inside/outside?
emmeans::emmeans(modelo_reg1, ~ Treatment|Planted) |>
multcomp::cld(adjust = "bonferroni")Planted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Strip_25 0.920 0.292 Inf 0.170 1.67 1
Nuclei_50 0.971 0.282 Inf 0.245 1.70 1
Full_Planted 1.162 0.249 Inf 0.520 1.80 1
Strip_50 1.278 0.261 Inf 0.606 1.95 1
Nuclei_25 1.280 0.280 Inf 0.558 2.00 1
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Nuclei_50 0.820 0.284 Inf 0.109 1.53 1
Strip_25 1.419 0.259 Inf 0.772 2.07 1
Strip_50 1.744 0.238 Inf 1.150 2.34 1
Nuclei_25 1.787 0.283 Inf 1.079 2.49 1
Full_Planted nonEst NA NA NA NA
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
P value adjustment: bonferroni method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Answer: No.
Is there any difference between inside/outside considering treatments?
emmeans::emmeans(modelo_reg1, ~ Planted|Treatment) |>
multcomp::cld(adjust = "bonferroni")Treatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.162 0.249 Inf 0.674 1.65 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.280 0.280 Inf 0.652 1.91 1
Outside 1.787 0.283 Inf 1.152 2.42 1
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 0.820 0.284 Inf 0.182 1.46 1
Inside 0.971 0.282 Inf 0.339 1.60 1
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 0.920 0.292 Inf 0.267 1.57 1
Outside 1.419 0.259 Inf 0.838 2.00 1
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.278 0.261 Inf 0.694 1.86 1
Outside 1.744 0.238 Inf 1.211 2.28 1
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Answer: No.
sjPlot::tab_model(modelo_reg1)Model matrix is rank deficient. Parameters
`TreatmentStrip_50:PlantedOutside` were not estimable.
| n | ||||
|---|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p | |
| Count Model | ||||
| (Intercept) | 3.20 | 1.96 – 5.21 | <0.001 | |
| Treatment [Nuclei_25] | 1.13 | 0.55 – 2.31 | 0.747 | |
| Treatment [Nuclei_50] | 0.83 | 0.40 – 1.70 | 0.605 | |
| Treatment [Strip_25] | 0.79 | 0.38 – 1.64 | 0.520 | |
| Treatment [Strip_50] | 1.12 | 0.56 – 2.25 | 0.743 | |
| Planted [Outside] | 1.59 | 0.81 – 3.14 | 0.178 | |
| Treatment [Nuclei_25] × Planted [Outside] |
1.04 | 0.38 – 2.82 | 0.937 | |
| Treatment [Nuclei_50] × Planted [Outside] |
0.54 | 0.20 – 1.49 | 0.233 | |
| Treatment [Strip_25] × Planted [Outside] |
1.03 | 0.38 – 2.80 | 0.949 | |
| (Intercept) | 3.53 | 2.03 – 9.52 | ||
| Zero-Inflated Model | ||||
| (Intercept) | 0.14 | 0.04 – 0.50 | 0.003 | |
| Random Effects | ||||
| σ2 | ||||
| τ00 Block | 0.00 | |||
| N Block | 5 | |||
| Observations | 180 | |||
broom.mixed::tidy(modelo_reg1) |>
dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |>
knitr::kable()| effect | component | group | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 1.162 | 0.249 | 4.663 | 0.000 |
| fixed | cond | NA | TreatmentNuclei_25 | 0.119 | 0.367 | 0.323 | 0.747 |
| fixed | cond | NA | TreatmentNuclei_50 | -0.191 | 0.369 | -0.517 | 0.605 |
| fixed | cond | NA | TreatmentStrip_25 | -0.241 | 0.376 | -0.643 | 0.520 |
| fixed | cond | NA | TreatmentStrip_50 | 0.116 | 0.355 | 0.328 | 0.743 |
| fixed | cond | NA | PlantedOutside | 0.466 | 0.346 | 1.347 | 0.178 |
| fixed | cond | NA | TreatmentNuclei_25:PlantedOutside | 0.040 | 0.508 | 0.080 | 0.937 |
| fixed | cond | NA | TreatmentNuclei_50:PlantedOutside | -0.617 | 0.518 | -1.192 | 0.233 |
| fixed | cond | NA | TreatmentStrip_25:PlantedOutside | 0.033 | 0.508 | 0.064 | 0.949 |
| fixed | cond | NA | TreatmentStrip_50:PlantedOutside | NA | NA | NA | NA |
| fixed | zi | NA | (Intercept) | -1.962 | 0.652 | -3.007 | 0.003 |
| ran_pars | cond | Block | sd__(Intercept) | 0.000 | NA | NA | NA |
Getting proportions of planted/unplanted species
prop_regenerant_part <- species_table |>
dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
dplyr::mutate(
planted_tree = ifelse(`Planting group`=="NA", "No", "Yes"))|>
dplyr::filter(n_sapling>0)|>
dplyr::group_by(Treatment, Planted, Block, Plot) |>
dplyr::summarise(
n_species = dplyr::n(),
n_ind = sum(n_sapling),
n_planted_ind = sum(ifelse(planted_tree == "Yes", n_sapling, 0), na.rm = TRUE),
n_planted_ind = ifelse(is.na(n_planted_ind), 0, n_planted_ind),
n_planted_spp = sum(planted_tree == "Yes"),
n_planted_spp = ifelse(is.na(n_planted_spp), 0, n_planted_spp),
.groups = "drop") |>
dplyr::mutate(
prop_planted_ind = n_planted_ind / n_ind,
prop_planted_spp = n_planted_spp / n_species
)
prop_regenerant_part |>
dplyr::group_by(Treatment, Planted)|>
dplyr::summarise(
prop_planted_mean = mean(prop_planted_ind, na.rm=T),
prop_planted_min = min(prop_planted_ind, na.rm=T),
prop_planted_max = max(prop_planted_ind, na.rm=T),
prop_planted_spp_mean = mean(prop_planted_spp, na.rm=T),
prop_planted_spp_min = min(prop_planted_spp, na.rm=T),
prop_planted_spp_max = max(prop_planted_spp, na.rm=T),
.groups = "drop"
) |>
dplyr::mutate(
prop_planted = sprintf("%.1f (%.1f - %.1f)",
round(prop_planted_mean,3)*100,
round(prop_planted_min,3)*100,
round(prop_planted_max,3)*100),
prop_planted_spp = sprintf("%.1f (%.1f - %.1f)",
round(prop_planted_spp_mean,3)*100,
round(prop_planted_spp_min,3)*100,
round(prop_planted_spp_max,3)*100)
) |>
dplyr::select(Treatment,Planted, prop_planted, prop_planted_spp) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = c(prop_planted, prop_planted_spp)
) # |> writexl::write_xlsx("resulta_prop_recruit.xlsx")# A tibble: 2 × 11
Planted prop_planted_Full_Plan…¹ prop_planted_Nuclei_25 prop_planted_Nuclei_50
<chr> <chr> <chr> <chr>
1 Inside 54.4 (0.0 - 100.0) 41.4 (0.0 - 100.0) 37.5 (0.0 - 100.0)
2 Outside <NA> 25.0 (0.0 - 91.2) 27.0 (0.0 - 92.3)
# ℹ abbreviated name: ¹prop_planted_Full_Planted
# ℹ 7 more variables: prop_planted_Strip_25 <chr>, prop_planted_Strip_50 <chr>,
# prop_planted_spp_Full_Planted <chr>, prop_planted_spp_Nuclei_25 <chr>,
# prop_planted_spp_Nuclei_50 <chr>, prop_planted_spp_Strip_25 <chr>,
# prop_planted_spp_Strip_50 <chr>
prop_regenerant_part |>
dplyr::group_by(Planted)|>
dplyr::summarise(
prop_planted_mean = mean(prop_planted_ind, na.rm=T),
prop_planted_se = (sd(prop_planted_ind, na.rm=T))/sqrt(dplyr::n()),
prop_planted_spp_mean = mean(prop_planted_spp, na.rm=T),
prop_planted_spp_se = (sd(prop_planted_spp, na.rm=T))/sqrt(dplyr::n()),
.groups = "drop"
)# A tibble: 2 × 5
Planted prop_planted_mean prop_planted_se prop_planted_spp_mean
<chr> <dbl> <dbl> <dbl>
1 Inside 0.524 0.0372 0.488
2 Outside 0.283 0.0309 0.243
# ℹ 1 more variable: prop_planted_spp_se <dbl>
Statistical analysis for proportion of individuals:
mod_binomial_part <-lme4::glmer(
cbind(n_planted_ind, n_ind - n_planted_ind) ~ Treatment*Planted + (1 | Block),
family = binomial(link = "logit"),
data = prop_regenerant_part)fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
performance::check_overdispersion(mod_binomial_part) # yes# Overdispersion test
dispersion ratio = 1.837
p-value = < 0.001
Overdispersion detected.
It should be Beta
mod_beta_part <-glmmTMB::glmmTMB(
cbind(n_planted_ind, n_ind - n_planted_ind) ~ Treatment*Planted + (1 | Block),
family = glmmTMB::betabinomial(link = "logit"),
data = prop_regenerant_part)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
performance::check_model(mod_beta_part)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
`check_outliers()` does not yet support models of class `glmmTMB`.
DHARMa::plotResiduals(mod_beta_part)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
Pairwise:
emmeans::emmeans(mod_beta_part, pairwise ~ Treatment | Planted) |> multcomp::cld(adjust="sidak") # no differencesPlanted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Nuclei_50 -0.429 0.311 Inf -1.228 0.370 1
Nuclei_25 -0.358 0.314 Inf -1.165 0.449 1
Full_Planted 0.243 0.326 Inf -0.595 1.082 1
Strip_25 0.447 0.298 Inf -0.319 1.212 1
Strip_50 0.490 0.292 Inf -0.260 1.240 1
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Strip_25 -1.196 0.298 Inf -1.938 -0.455 1
Nuclei_25 -1.096 0.287 Inf -1.810 -0.382 1
Nuclei_50 -0.935 0.314 Inf -1.718 -0.153 1
Strip_50 -0.376 0.264 Inf -1.035 0.282 1
Full_Planted nonEst NA NA NA NA
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: sidak method for varying numbers of estimates
Results are given on the log odds ratio (not the response) scale.
P value adjustment: sidak method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
emmeans::emmeans(mod_beta_part, pairwise ~ Planted | Treatment) |> multcomp::cld(adjust="sidak") # All differentTreatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 0.243 0.326 Inf -0.397 0.883 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -1.096 0.287 Inf -1.737 -0.455 1
Inside -0.358 0.314 Inf -1.061 0.345 1
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -0.935 0.314 Inf -1.638 -0.233 1
Inside -0.429 0.311 Inf -1.125 0.267 1
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -1.196 0.298 Inf -1.862 -0.530 1
Inside 0.447 0.298 Inf -0.220 1.113 2
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -0.376 0.264 Inf -0.968 0.215 1
Inside 0.490 0.292 Inf -0.163 1.143 2
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: sidak method for varying numbers of estimates
Results are given on the log odds ratio (not the response) scale.
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Statistical analysis for proportion of individuals:
mod_binomial_prop <-lme4::glmer(
cbind(n_planted_spp, n_species - n_planted_spp) ~ Treatment*Planted + (1 | Block),
family = binomial(link = "logit"),
data = prop_regenerant_part)fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
boundary (singular) fit: see help('isSingular')
performance::check_overdispersion(mod_binomial_prop) # underdispersion# Overdispersion test
dispersion ratio = 0.769
p-value = 0.032
Underdispersion detected.
It should be Beta
mod_beta_prop <-glmmTMB::glmmTMB(
cbind(n_planted_spp, n_species - n_planted_spp) ~ Treatment*Planted + (1 | Block),
family = glmmTMB::betabinomial(link = "logit"),
data = prop_regenerant_part)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
problem; false convergence (8). See vignette('troubleshooting'),
help('diagnose')
performance::check_model(mod_beta_prop)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
`check_outliers()` does not yet support models of class `glmmTMB`.
DHARMa::plotResiduals(mod_beta_prop)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
Pairwise:
emmeans::emmeans(mod_beta_prop, pairwise ~ Treatment | Planted) |> multcomp::cld(adjust="sidak") # no differencesPlanted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Nuclei_25 -0.582 0.286 Inf -1.318 0.154 1
Nuclei_50 -0.528 0.263 Inf -1.203 0.147 1
Strip_25 -0.210 0.246 Inf -0.841 0.421 1
Full_Planted 0.251 0.291 Inf -0.496 0.999 1
Strip_50 0.297 0.259 Inf -0.368 0.962 1
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Nuclei_25 -1.399 0.250 Inf -2.021 -0.777 1
Strip_25 -1.204 0.233 Inf -1.784 -0.624 1
Nuclei_50 -0.944 0.257 Inf -1.585 -0.304 1
Strip_50 -0.847 0.208 Inf -1.366 -0.329 1
Full_Planted nonEst NA NA NA NA
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: sidak method for varying numbers of estimates
Results are given on the log odds ratio (not the response) scale.
P value adjustment: sidak method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
emmeans::emmeans(mod_beta_prop, pairwise ~ Planted | Treatment) |> multcomp::cld(adjust="sidak") # All differentTreatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 0.251 0.291 Inf -0.319 0.8216 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -1.399 0.250 Inf -1.957 -0.8403 1
Inside -0.582 0.286 Inf -1.222 0.0587 2
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -0.944 0.257 Inf -1.520 -0.3693 1
Inside -0.528 0.263 Inf -1.116 0.0599 1
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -1.204 0.233 Inf -1.724 -0.6835 1
Inside -0.210 0.246 Inf -0.759 0.3398 2
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -0.847 0.208 Inf -1.313 -0.3820 1
Inside 0.297 0.259 Inf -0.282 0.8763 2
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: sidak method for varying numbers of estimates
Results are given on the log odds ratio (not the response) scale.
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
broom.mixed::tidy(mod_beta_part) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 11 × 8
effect component group term estimate std.error statistic p.value
<chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 fixed cond <NA> (Intercept) 0.243 0.326 0.745 0.456
2 fixed cond <NA> TreatmentNucle… -0.601 0.454 -1.32 0.185
3 fixed cond <NA> TreatmentNucle… -0.672 0.451 -1.49 0.136
4 fixed cond <NA> TreatmentStrip… 0.203 0.442 0.461 0.645
5 fixed cond <NA> TreatmentStrip… 0.247 0.437 0.564 0.572
6 fixed cond <NA> PlantedOutside -0.867 0.394 -2.20 0.028
7 fixed cond <NA> TreatmentNucle… 0.129 0.577 0.223 0.824
8 fixed cond <NA> TreatmentNucle… 0.36 0.592 0.608 0.543
9 fixed cond <NA> TreatmentStrip… -0.776 0.576 -1.35 0.178
10 fixed cond <NA> TreatmentStrip… NA NA NA NA
11 ran_pars cond Block sd__(Intercept) 0 NA NA NA
sjPlot::tab_model(mod_beta_part)Model matrix is rank deficient. Parameters
`TreatmentStrip_50:PlantedOutside` were not estimable.
| cbind(n planted ind,n ind-n planted ind) |
|||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 1.28 | 0.67 – 2.42 | 0.456 |
| Treatment [Nuclei_25] | 0.55 | 0.23 – 1.33 | 0.185 |
| Treatment [Nuclei_50] | 0.51 | 0.21 – 1.24 | 0.136 |
| Treatment [Strip_25] | 1.23 | 0.52 – 2.91 | 0.645 |
| Treatment [Strip_50] | 1.28 | 0.54 – 3.02 | 0.572 |
| Planted [Outside] | 0.42 | 0.19 – 0.91 | 0.028 |
| Treatment [Nuclei_25] × Planted [Outside] |
1.14 | 0.37 – 3.52 | 0.824 |
| Treatment [Nuclei_50] × Planted [Outside] |
1.43 | 0.45 – 4.57 | 0.543 |
| Treatment [Strip_25] × Planted [Outside] |
0.46 | 0.15 – 1.42 | 0.178 |
| Random Effects | |||
| σ2 | 0.37 | ||
| τ00 Block | 0.00 | ||
| N Block | 5 | ||
| Observations | 176 | ||
| Marginal R2 / Conditional R2 | 0.500 / NA | ||
Estimating density of recruits (\(ind.ha^{-1}\)) and without ruderal (Baccharis and Vernonanthura)
density_ruderal_n_m2 <- species_table |>
dplyr::filter(scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
# dplyr::select(Treatment, Block, Plot, Planted)|>unique()|> with(table(Block, Plot, by=Treatment)) # Full_Planted P2B5 & P3B1 (Inside)
dplyr::add_row(
scientificName = c(NA, NA),
`Seed dispersal` = c(NA, NA),
`Planting group` = c(NA, NA),
Treatment = c("Full_Planted", "Full_Planted"),
Block = c("B5", "B1"),
Plot = c("P2", "P3"),
Planted = c("Inside","Inside"),
n_sapling= c(0,0),
n_tree = c(0,0)
) |>
dplyr::mutate(
n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
dplyr::group_by(Treatment, Block, Planted) |>
dplyr::summarise(
n_small_m2 = (sum(n_small, na.rm=T))/(144*4),
n_large_m2 = (sum(n_large, na.rm=T))/(144*4),
.groups = "drop")
dens_ruderal_ha <- density_ruderal_n_m2 |>
dplyr::mutate(
n_small_ha = n_small_m2 * 10000,
n_large_ha = n_large_m2 * 10000
)
dens_ruderal_P_UP <- dens_ruderal_ha |>
tidyr::pivot_wider(
names_from = Planted,
values_from = c(n_small_ha, n_large_ha),
names_sep = "_"
)
prop_table <- tidyr::tibble(
Treatment = unique(sort(species_table$Treatment)),
prop_P = c(1, 0.25, 0.5, 0.25, 0.5),
prop_UP = 1 - prop_P
)
result_ruderal <- dens_ruderal_P_UP |>
dplyr::left_join(prop_table, by = "Treatment") |>
dplyr::mutate(
n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0) + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
) |>
dplyr::select(Treatment,
starts_with("n_small_ha"),
starts_with("n_large_ha"),
n_small_total_ha,
n_large_total_ha
)Smaller recruits
result_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_ha_Inside, na.rm=TRUE),1),
min = round(min(n_small_ha_Inside, na.rm=TRUE),0),
max = round(max(n_small_ha_Inside, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 2545.1 (1580.0 - 4566.0) 4388.9 (2847.0 - 6858.0) 5145.8 (1… 4093.8 … 5409.7 …
result_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
min = round(min(n_small_ha_Outside, na.rm=TRUE), 0),
max = round(max(n_small_ha_Outside, na.rm=TRUE), 0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_small_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 NaN (Inf - -Inf) 2784.7 (1181.0 - 6858.0) 2638.9 (1302.0 - … 2246.5 … 2489.6 …
result_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_total_ha, na.rm=TRUE),1),
min = round(min(n_small_total_ha, na.rm=TRUE),0),
max = round(max(n_small_total_ha, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 2545.1 (1580.0 - 4566.0) 1592.9 (712.0 - 5143.0) 1946.2 (65… 1354.2 … 1974.8 …
Larger recruits
All treatments and inside/outside areas showed similar values
result_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_ha_Inside, na.rm=TRUE),1),
min = round(min(n_large_ha_Inside, na.rm=TRUE),0),
max = round(max(n_large_ha_Inside, na.rm=TRUE),0)
) |>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 13.9 (0.0 - 35.0) 17.4 (0.0 - 69.0) 3.5 (0.0 - 17.0) 3.5 (0.0 - 17.0) 3.5 (0.…
result_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_ha_Outside, na.rm=TRUE),1),
min = round(min(n_large_ha_Outside, na.rm=TRUE),0),
max = round(max(n_large_ha_Outside, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_large_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 NaN (Inf - -Inf) 31.2 (0.0 - 52.0) 24.3 (0.0 - 87.0) 48.6 (0.0 - 139… 17.4 (0…
result_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_total_ha, na.rm = TRUE), 1),
min = round(min(n_large_total_ha, na.rm = TRUE), 0),
max = round(max(n_large_total_ha, na.rm = TRUE), 0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 13.9 (0.0 - 35.0) 13.9 (0.0 - 39.0) 6.9 (0.0 - 43.0) 18.7 (0.0 - 104… 5.2 (0.…
ruderal_df <- species_table |>
dplyr::filter(scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
# dplyr::select(Treatment, Block, Plot, Planted)|>unique()|> with(table(Block, Plot, by=Treatment)) # Full_Planted P2B5 & P3B1 (Inside)
dplyr::add_row(
scientificName = c(NA, NA),
`Seed dispersal` = c(NA, NA),
`Planting group` = c(NA, NA),
Treatment = c("Full_Planted", "Full_Planted"),
Block = c("B5", "B1"),
Plot = c("P2", "P3"),
Planted = c("Inside","Inside"),
n_sapling= c(0,0),
n_tree = c(0,0),
n_total= c(0,0)
)|>
dplyr::group_by(Treatment, Block, Plot, Planted)|>
dplyr::summarise(n = sum(n_total), .groups="drop")Plotting
ruderal_df |>
ggplot(aes(x=Planted, y=n, color=Block))+
geom_point()+facet_wrap(~Treatment)+theme_bw()Models
modelo_rud1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = ruderal_df)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
modelo_rud2 <- glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = ruderal_df)
modelo_rud3 <- glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = ruderal_df)
modelo_rud_null <- glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
ziformula = ~1,
# family = poisson(), # ZIP
family = glmmTMB::nbinom2(), # ZINB
# family = glmmTMB::truncated_nbinom2(), # hurdle
data = ruderal_df)
# Model (null) | (treat*planted) | (treat) | (planted)
# ZIP
AIC(modelo_rud_null, modelo_rud1, modelo_rud2, modelo_rud3) df AIC
modelo_rud_null 3 1725.166
modelo_rud1 12 1658.610
modelo_rud2 8 1707.148
modelo_rud3 5 1660.457
performance::check_overdispersion(modelo_rud1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
# Overdispersion test
dispersion ratio = 0.965
p-value = 0.872
No overdispersion detected.
performance::check_model(modelo_rud1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
`check_outliers()` does not yet support models of class `glmmTMB`.
Exploring residuals
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_rud1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
DHARMa::plotQQunif(modelo_rud1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
DHARMa::plotResiduals(modelo_rud3)
DHARMa::plotQQunif(modelo_rud3)car::Anova(modelo_rud1)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: n
Chisq Df Pr(>Chisq)
Treatment 16.1644 4 0.002806 **
Planted 65.9720 1 4.574e-16 ***
Treatment:Planted 1.5434 3 0.672287
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Is there any difference between treatments considering inside/outside?
emmeans::emmeans(modelo_rud1, ~ Treatment|Planted) |>
multcomp::cld(adjust = "bonferroni")Planted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Full_Planted 3.66 0.206 Inf 3.13 4.19 1
Strip_25 4.09 0.201 Inf 3.57 4.60 12
Nuclei_25 4.13 0.201 Inf 3.61 4.65 12
Nuclei_50 4.23 0.200 Inf 3.72 4.75 2
Strip_50 4.29 0.200 Inf 3.78 4.81 2
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Strip_25 3.42 0.203 Inf 2.91 3.93 1
Strip_50 3.48 0.202 Inf 2.97 3.99 1
Nuclei_50 3.51 0.202 Inf 3.00 4.01 1
Nuclei_25 3.60 0.204 Inf 3.09 4.11 1
Full_Planted nonEst NA NA NA NA
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
P value adjustment: bonferroni method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Answer: Yes (planted) and no (only unplanted areas)
Is there any difference between inside/outside considering treatments?
emmeans::emmeans(modelo_rud1, ~ Planted|Treatment) |>
multcomp::cld(adjust = "bonferroni")Treatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 3.66 0.206 Inf 3.26 4.06 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.60 0.204 Inf 3.15 4.06 1
Inside 4.13 0.201 Inf 3.68 4.58 2
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.51 0.202 Inf 3.05 3.96 1
Inside 4.23 0.200 Inf 3.79 4.68 2
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.42 0.203 Inf 2.97 3.88 1
Inside 4.09 0.201 Inf 3.64 4.54 2
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.48 0.202 Inf 3.03 3.93 1
Inside 4.29 0.200 Inf 3.85 4.74 2
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Answer: No.
sjPlot::tab_model(modelo_rud1)Model matrix is rank deficient. Parameters
`TreatmentStrip_50:PlantedOutside` were not estimable.
| n | ||||
|---|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p | |
| Count Model | ||||
| (Intercept) | 38.84 | 25.95 – 58.15 | <0.001 | |
| Treatment [Nuclei_25] | 1.60 | 1.14 – 2.24 | 0.006 | |
| Treatment [Nuclei_50] | 1.78 | 1.27 – 2.49 | 0.001 | |
| Treatment [Strip_25] | 1.53 | 1.09 – 2.15 | 0.013 | |
| Treatment [Strip_50] | 1.89 | 1.35 – 2.64 | <0.001 | |
| Planted [Outside] | 0.44 | 0.32 – 0.61 | <0.001 | |
| Treatment [Nuclei_25] × Planted [Outside] |
1.33 | 0.84 – 2.12 | 0.228 | |
| Treatment [Nuclei_50] × Planted [Outside] |
1.09 | 0.69 – 1.73 | 0.711 | |
| Treatment [Strip_25] × Planted [Outside] |
1.16 | 0.73 – 1.85 | 0.529 | |
| (Intercept) | 51.38 | 23.07 – 140.38 | ||
| Zero-Inflated Model | ||||
| (Intercept) | 0.01 | 0.00 – 0.05 | <0.001 | |
| Random Effects | ||||
| σ2 | 0.02 | |||
| τ00 Block | 0.13 | |||
| ICC | 0.87 | |||
| N Block | 5 | |||
| Observations | 179 | |||
| Marginal R2 / Conditional R2 | 0.424 / 0.924 | |||
broom.mixed::tidy(modelo_rud1) |>
dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |>
dplyr::select(-c(component, group))|>
knitr::kable()| effect | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|
| fixed | (Intercept) | 3.659 | 0.206 | 17.777 | 0.000 |
| fixed | TreatmentNuclei_25 | 0.471 | 0.171 | 2.748 | 0.006 |
| fixed | TreatmentNuclei_50 | 0.575 | 0.171 | 3.356 | 0.001 |
| fixed | TreatmentStrip_25 | 0.427 | 0.172 | 2.486 | 0.013 |
| fixed | TreatmentStrip_50 | 0.635 | 0.171 | 3.713 | 0.000 |
| fixed | PlantedOutside | -0.814 | 0.167 | -4.873 | 0.000 |
| fixed | TreatmentNuclei_25:PlantedOutside | 0.286 | 0.237 | 1.207 | 0.228 |
| fixed | TreatmentNuclei_50:PlantedOutside | 0.087 | 0.236 | 0.371 | 0.711 |
| fixed | TreatmentStrip_25:PlantedOutside | 0.149 | 0.237 | 0.630 | 0.529 |
| fixed | TreatmentStrip_50:PlantedOutside | NA | NA | NA | NA |
| fixed | (Intercept) | -4.497 | 0.721 | -6.238 | 0.000 |
| ran_pars | sd__(Intercept) | 0.365 | NA | NA | NA |
Small letters: Comparison between treatments, only inside planted areas.
Capital letters: Comparison between treatments, only outside planted areas.
ordem <- c("Nuclei_25", "Strip_25", "Nuclei_50", "Strip_50", "Full_Planted")
recruits_ha <- species_table |>
dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes")) |>
dplyr::mutate(
n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
dplyr::group_by(Treatment, Block, Planted) |>
dplyr::summarise(
n_small_m2 = (sum(n_small, na.rm=T))/(144*4),
n_large_m2 = (sum(n_large, na.rm=T))/(144*4),
.groups = "drop") |>
dplyr::mutate(
n_small_ha = n_small_m2 * 10000,
n_large_ha = n_large_m2 * 10000
) |>
tidyr::pivot_wider(
names_from = Planted,
values_from = c(n_small_ha, n_large_ha),
names_sep = "_"
) |>
dplyr::select(Treatment, Block, n_small_ha_Inside:n_large_ha_Outside)
ruderal_ha <- species_table |>
dplyr::filter(scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes")) |>
dplyr::mutate(
n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
dplyr::group_by(Treatment, Block, Planted) |>
dplyr::summarise(
n_small_m2 = (sum(n_small, na.rm=T))/(144*4),
n_large_m2 = (sum(n_large, na.rm=T))/(144*4),
.groups = "drop") |>
dplyr::mutate(
n_small_ha = n_small_m2 * 10000,
n_large_ha = n_large_m2 * 10000
) |>
tidyr::pivot_wider(
names_from = Planted,
values_from = c(n_small_ha, n_large_ha),
names_sep = "_"
) |>
dplyr::select(Treatment, Block, n_small_ha_Inside:n_large_ha_Outside)
prop_table <- tidyr::tibble(
Treatment = unique(sort(species_table$Treatment)),
prop_P = c(1, 0.25, 0.5, 0.25, 0.5),
prop_UP = 1 - prop_P
)
data_fig_recruits <- recruits_ha |>
dplyr::left_join(prop_table, by = "Treatment") |>
dplyr::mutate(
n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0) + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
) |>
dplyr::select(Treatment,
starts_with("n_small_ha"),
starts_with("n_large_ha"),
n_small_total_ha,
n_large_total_ha
)
data_fig_ruderal <- ruderal_ha |>
dplyr::left_join(prop_table, by = "Treatment") |>
dplyr::mutate(
n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0) + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
) |>
dplyr::select(Treatment,
starts_with("n_small_ha"),
starts_with("n_large_ha"),
n_small_total_ha,
n_large_total_ha
)Small regenerants
data_fig_small_1 <- data_fig_recruits |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean_inside = round(mean(n_small_ha_Inside, na.rm=TRUE), 1),
se_inside = round((sd(n_small_ha_Inside, na.rm=TRUE))/sqrt(dplyr::n()), 1),
mean_outside = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
se_outside = round((sd(n_small_ha_Outside, na.rm=TRUE))/sqrt(dplyr::n()), 1)
) |>
tidyr::pivot_longer(cols = -Treatment, names_to = c(".value", "side"), names_sep = "_") |>
dplyr::filter(!is.na(se)) |>
dplyr::mutate(class="regenerant")
data_fig_small_2 <- data_fig_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean_inside = round(mean(n_small_ha_Inside, na.rm=TRUE), 1),
se_inside = round((sd(n_small_ha_Inside, na.rm=TRUE))/sqrt(dplyr::n()), 1),
mean_outside = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
se_outside = round((sd(n_small_ha_Outside, na.rm=TRUE))/sqrt(dplyr::n()), 1)
) |>
tidyr::pivot_longer(cols = -Treatment, names_to = c(".value", "side"), names_sep = "_") |>
dplyr::filter(!is.na(se)) |>
dplyr::mutate(class="ruderal")
data_fig_small <- dplyr::bind_rows(data_fig_small_1,data_fig_small_2)New Figure
ordem_x <- c(
"Nuclei_25_inside", "Nuclei_25_outside",
"Strip_25_inside", "Strip_25_outside",
"Nuclei_50_inside", "Nuclei_50_outside",
"Strip_50_inside", "Strip_50_outside",
"Full_Planted_inside"
)
data_fig_small_top2 <- data_fig_small |>
group_by(Treatment, side) |>
arrange(desc(class)) |>
mutate(cum_mean = cumsum(mean)) |>
ungroup() |>
mutate(Treat_side = interaction(Treatment, side, sep = "_"),
Treat_side = factor(Treat_side, levels = ordem_x))
letras <- data.frame(
Treat_side = levels(data_fig_small_top2$Treat_side),
label = c("Ab","Aa",
"Ab","Aa",
"Aa","Aa",
"Ab","Aa",
"A"),
y = c(5000, 4500,
5000, 4000,
6500, 4000,
6500, 4500,
3500)
)
fig_small<-ggplot(data_fig_small_top2,
aes(x = Treat_side, y = mean, fill = side, pattern = class)) +
geom_col_pattern(
position = "stack",
color = "black",
pattern_fill = "#92D050", pattern_alpha = 0.5,
pattern_angle = 45, pattern_density = 0.1,
pattern_spacing = 0.1
) +
scale_pattern_manual(
values = c("regenerant" = "none", "ruderal" = "stripe"),
guide = "none") +
scale_fill_manual(labels = c("Planted", "Unplanted"),
values = c("#92D050", "white")) +
geom_errorbar(
data = data_fig_small_top2 |> group_by(Treatment, side) |> slice_tail(n = 1),
aes(ymin = cum_mean - se, ymax = cum_mean + se, x = Treat_side),
width = .2,
inherit.aes = FALSE
) +
labs(fill = "",
y = expression("Recruits (n.ha"^{-1}~")"),
x = "") +
geom_text(data = letras,
aes(x = Treat_side, y = y, label = label),
inherit.aes = FALSE,
size = 3.5, vjust = -0.5) +
scale_x_discrete(labels = function(x) gsub("_", "\n", x)) +
scale_y_continuous(expand =c(0,0), limits = c(0,7500), breaks = seq(0, 7000, by=1500))+
theme_classic() +
theme(
axis.text.x = element_blank(), #element_text(color="black", hjust = 0.5, size=9),
axis.text.y = element_text(color="black", size=10),
axis.title.y = element_text(color="black", size=11),
legend.position = c(0.18, 0.95),
legend.text = element_text(size = 10),
legend.key.size = unit(3, "mm"),
legend.spacing = unit(0, "mm"),
legend.background = element_blank(),
legend.title = element_text(size = 8))
fig_small# ggsave("fig_small_recruits.png", plot=fig_small, width=90, height = 90 ,unit="mm")Without overabundant species
data_fig_recruits |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean_inside = round(mean(n_small_ha_Inside, na.rm=TRUE), 1),
se_inside = round((sd(n_small_ha_Inside, na.rm=TRUE))/sqrt(dplyr::n()), 1),
mean_outside = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
se_outside = round((sd(n_small_ha_Outside, na.rm=TRUE))/sqrt(dplyr::n()), 1)
) |>
tidyr::pivot_longer(cols = -Treatment, names_to = c(".value", "side"), names_sep = "_") |>
dplyr::filter(!is.na(se)) |>
dplyr::mutate(
tratamento = factor(Treatment, ordered = TRUE, levels = c("Nuclei_25", "Strip_25", "Nuclei_50", "Strip_50","Full_Planted")),
tratamento = forcats::fct_recode(tratamento,
"Nuclei\n25%" = "Nuclei_25", "Strips\n25%" = "Strip_25",
"Nuclei\n50%" = "Nuclei_50","Strips\n50%" = "Strip_50","Planting\n100%" = "Full_Planted")
)|>
ggplot(aes(x = tratamento, y = mean, fill = side)) +
geom_bar(position = position_dodge2(preserve = "single"),
stat = "identity", color = "black") +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = .2, position = position_dodge(0.9)) +
scale_fill_manual(labels = c('Planted', 'Unplanted'),
values = c("#92D050", "white")) +
# geom_text(data = letras, aes(x = tratamento, y = y, label = label), size=3.5, vjust = -0.5) +
scale_y_continuous(expand =c(0,0), limits = c(0,1700), breaks = seq(0, 1700, by=500))+
labs(fill = "",
y = expression("Recruits (n.ha"^{-1}~")"),
x = "") +
theme_classic()+
theme(
axis.text.x = element_text(color="black", hjust = 0.5, size=10),
axis.text.y = element_text(color="black", size=10),
axis.title.y = element_text(color="black", size=11),
legend.position = c(0.16, 0.98),
plot.margin = margin(t=5, l=2, r=2, 0, "mm"),
legend.text = element_text(size = 10, margin = margin(t=0,0,0,l=1, "mm")),
legend.key.size = unit(3, "mm"),
legend.spacing = unit(0, "mm"),
legend.background = element_blank(),
legend.title = element_text(size = 8))data_fig_large <- result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean_inside = round(mean(n_large_ha_Inside, na.rm=TRUE), 1),
se_inside = round((sd(n_large_ha_Inside, na.rm=TRUE))/sqrt(dplyr::n()), 1),
mean_outside = round(mean(n_large_ha_Outside, na.rm=TRUE), 1),
se_outside = round((sd(n_large_ha_Outside, na.rm=TRUE))/sqrt(dplyr::n()), 1)
) |>
tidyr::pivot_longer(cols = -Treatment, names_to = c(".value", "side"), names_sep = "_") |>
dplyr::mutate(
Treatment = ordered(Treatment,levels = ordem),
tratamento = forcats::fct_recode(Treatment,
"Nuclei 25%" = "Nuclei_25", "Strips 25%" = "Strip_25",
"Nuclei 50%" = "Nuclei_50", "Strips 50%" = "Strip_50","Planting 100%" = "Full_Planted"),
Treat_side = interaction(Treatment, side, sep = "_"),
Treat_side = factor(Treat_side, ordered = TRUE, levels = ordem_x,
labels = c("Nuclei\n 25%", " " ,
"Strips\n 25%", " ",
"Nuclei\n 50%", " ",
"Strips\n 50%", " ",
"Planting\n100%"))
)|>
dplyr::filter(!is.na(se))Inside and outside planted areas are similar
letras <- data.frame(
tratamento = rep(levels(as.factor(data_fig_large$tratamento)), each=2),
side = levels(as.factor(data_fig_large$side)),
label = c("A "," a",
"A "," a",
"A "," a",
"A "," a",
"A",""),
y = c(300, 470,
200, 340,
200, 200,
300, 470,
320, 400)
)
fig_large <- data_fig_large |>
ggplot(aes(x = tratamento, y = mean, fill = side)) +
geom_bar(
position = position_dodge2(width = 1, padding = .1,preserve = "single"),
stat = "identity", color = "black") +
geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = .2, position = position_dodge(0.9)) +
scale_fill_manual(labels = c('Planted', 'Unplanted'), values = c("#92D050", "white")) +
geom_text(data = letras, aes(x = tratamento, y = y, label = label), size=3.5, vjust = -0.5) +
scale_y_continuous(expand =c(0,0), limits = c(0,550), breaks = seq(0, 500, by=100))+
labs(fill = "",
y = expression("Recruits (n.ha"^{-1}*")"),
x = "") +
theme_classic()+
theme(
axis.text.x = element_text(color="black", hjust = 0.5, size=10),
axis.text.y = element_text(color="black", size=10),
axis.title.y = element_text(color="black", size=11),
plot.margin = margin(t=5, l=2, r=2,0, "mm"),
legend.position = "none", # c(0.90, 0.95),
legend.text = element_text(size = 9, margin = margin(0,0,0,l=1, "mm")),
legend.key.size = unit(3, "mm"),
legend.spacing = unit(0, "mm"),
legend.background = element_blank(),
legend.title = element_text(size = 8))
fig_largeAnother way
letras <- data.frame(
Treat_side = levels(data_fig_large$Treat_side),
label = c("Aa","Aa",
"Aa","Aa",
"Aa","Aa",
"Aa","Aa",
"Aa"),
y = c(300, 470,
200, 340,
200, 200,
300, 470,
320)
)
data_fig_large <- data_fig_large |>
dplyr::mutate(hjust_value = ifelse(Treat_side == last(Treat_side), 0.5, 0))
fig_large<-ggplot(data_fig_large,
aes(x = Treat_side, y = mean, fill = side)) +
geom_bar(
stat = "identity",
position = "stack",
color = "black") +
scale_fill_manual(labels = c("Planted", "Unplanted"),
values = c("#92D050", "white")) +
geom_errorbar(
aes(ymin = mean - se, ymax = mean + se, x = Treat_side),
width = .2
) +
labs(fill = "",
y = expression("Recruits (n.ha"^{-1}*")"),
x = "") +
geom_text(data = letras,
aes(x = Treat_side, y = y, label = label),
inherit.aes = FALSE,
size = 3.5, vjust = -0.5) +
scale_y_continuous(expand =c(0,0), limits = c(0,550), breaks = seq(0, 500, by=100))+
theme_classic() +
theme(
axis.text.x = element_text(color="black", hjust = data_fig_large$hjust_value, size=9),
axis.text.y = element_text(color="black", size=10),
axis.title.y = element_text(color="black", size=11),
legend.position = "none",# c(0.18, 0.95),
legend.text = element_text(size = 10),
legend.key.size = unit(3, "mm"),
legend.spacing = unit(0, "mm"),
legend.background = element_blank(),
legend.title = element_text(size = 8))Warning: Vectorized input to `element_text()` is not officially supported.
ℹ Results may be unexpected or may change in future versions of ggplot2.
fig_largefig_recruits <- cowplot::plot_grid(fig_small + theme(axis.text.x=element_blank()),
fig_large,
ncol = 1,
labels = c("(a)", "(b)"), label_size = 10, label_x = 0.01, label_y = c(1,1.1)
)
fig_recruitsggsave("fig_recruits.png", plot=fig_recruits, width=90, height = 100 ,unit="mm")Only recruits (not including planted individuals)
table_recruits <- readxl::read_xlsx("Itatinga_Nuc_recruitment_SpeciesList.xlsx")
table_recruits |> knitr::kable(format = "html")| family | scientificName | Nuclei 25% (small) | Nuclei 25% (large) | Nuclei 25% (total) | Strips 25% (small) | Strips 25% (large) | Strips 25% (total) | Nuclei 50% (small) | Nuclei 50% (large) | Nuclei 50% (total) | Strips 50% (small) | Strips 50% (large) | Strips 50% (total) | Planting 100% (small) | Planting 100% (large) | Planting 100% (total) | Total small | Total large | Total |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Anacardiaceae | Astronium graveolens | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 5 | 5 | 0 | 0 | 0 | 1 | 5 | 6 |
| Anacardiaceae | Myracrodruon urundeuva | 0 | 12 | 12 | 0 | 2 | 2 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 16 | 16 |
| Anacardiaceae | Schinus terebinthifolia | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 |
| Anacardiaceae | Tapirira guianensis | 0 | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3 |
| Annonaceae | Annona cacans | 0 | 0 | 0 | 11 | 0 | 11 | 0 | 0 | 0 | 0 | 2 | 2 | 1 | 0 | 1 | 12 | 2 | 14 |
| Annonaceae | Annona coriacea | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 2 | 4 |
| Annonaceae | Duguetia furfuracea | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 1 | 2 | 1 | 3 |
| Apocynaceae | Aspidosperma tomentosum | 2 | 1 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 4 | 4 | 2 | 5 | 7 |
| Araliaceae | Didymopanax vinosus | 0 | 0 | 0 | 4 | 3 | 7 | 2 | 0 | 2 | 3 | 1 | 4 | 0 | 0 | 0 | 9 | 4 | 13 |
| Arecaceae | Syagrus romanzoffiana | 1 | 1 | 2 | 1 | 0 | 1 | 3 | 0 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 5 | 1 | 6 |
| Asteraceae | Baccharis dracunculifolia | 1432 | 14 | 1446 | 916 | 15 | 931 | 931 | 8 | 939 | 1784 | 6 | 1790 | 511 | 4 | 515 | 5574 | 47 | 5621 |
| Asteraceae | Moquiniastrum barrosoae | 0 | 0 | 0 | 3 | 1 | 4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 1 | 4 |
| Asteraceae | Moquiniastrum polymorphum | 0 | 17 | 17 | 2 | 3 | 5 | 13 | 2 | 15 | 2 | 0 | 2 | 0 | 0 | 0 | 17 | 22 | 39 |
| Asteraceae | Piptocarpha rotundifolia | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Asteraceae | Vernonanthura polyanthes | 634 | 0 | 634 | 910 | 0 | 910 | 1311 | 0 | 1311 | 491 | 0 | 491 | 0 | 2 | 2 | 3346 | 2 | 3348 |
| Bignoniaceae | Cybistax antisyphilitica | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| Bignoniaceae | Handroanthus albus | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 9 | 0 | 9 | 9 | 2 | 11 |
| Bignoniaceae | Handroanthus ochraceus | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 2 | 3 |
| Boraginaceae | Cordia sellowiana | 0 | 0 | 0 | 0 | 6 | 6 | 0 | 3 | 3 | 0 | 3 | 3 | 0 | 0 | 0 | 0 | 12 | 12 |
| Burseraceae | Protium spruceanum | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 1 |
| Calophyllaceae | Myrsine gardneriana | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Cannabaceae | Trema micrantha | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 10 | 10 | 0 | 0 | 0 | 0 | 11 | 11 |
| Chrysobalanaceae | Couepia grandiflora | 1 | 0 | 1 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 2 | 3 |
| Combretaceae | Terminalia argentea | 0 | 0 | 0 | 0 | 2 | 2 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 4 | 4 |
| Euphorbiaceae | Alchornea triplinervia | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 44 | 0 | 44 | 45 | 0 | 45 |
| Euphorbiaceae | Croton floribundus | 73 | 0 | 73 | 55 | 1 | 56 | 119 | 0 | 119 | 147 | 0 | 147 | 0 | 0 | 0 | 394 | 1 | 395 |
| Euphorbiaceae | Croton urucurana | 1 | 0 | 1 | 2 | 0 | 2 | 0 | 0 | 0 | 7 | 1 | 8 | 0 | 0 | 0 | 10 | 1 | 11 |
| Euphorbiaceae | Mabea fistulifera | 1 | 0 | 1 | 12 | 0 | 12 | 7 | 0 | 7 | 31 | 2 | 33 | 0 | 3 | 3 | 51 | 5 | 56 |
| Euphorbiaceae | Sapium glandulosum | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 2 | 0 | 0 | 0 | 1 | 2 | 3 |
| Fabaceae | Anadenanthera colubrina var. cebil | 2 | 5 | 7 | 5 | 0 | 5 | 0 | 1 | 1 | 8 | 13 | 21 | 0 | 0 | 0 | 15 | 19 | 34 |
| Fabaceae | Andira fraxinifolia | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Fabaceae | Bauhinia rufa | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 0 | 0 | 2 | 0 | 2 |
| Fabaceae | Chamaecrista comosa | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| Fabaceae | Copaifera langsdorffii | 2 | 2 | 4 | 2 | 4 | 6 | 2 | 2 | 4 | 5 | 0 | 5 | 0 | 0 | 0 | 11 | 8 | 19 |
| Fabaceae | Dalbergia frutescens | 0 | 2 | 2 | 7 | 0 | 7 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 7 | 2 | 9 |
| Fabaceae | Dalbergia miscolobium | 0 | 0 | 0 | 1 | 2 | 3 | 1 | 1 | 2 | 2 | 1 | 3 | 0 | 0 | 0 | 4 | 4 | 8 |
| Fabaceae | Enterolobium contortisiliquum | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Fabaceae | Erythrina mulungu | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 1 |
| Fabaceae | Hymenaea stigonocarpa | 0 | 2 | 2 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3 |
| Fabaceae | Inga vera | 0 | 0 | 0 | 4 | 0 | 4 | 4 | 0 | 4 | 3 | 0 | 3 | 1 | 0 | 1 | 12 | 0 | 12 |
| Fabaceae | Leptolobium dasycarpum | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3 |
| Fabaceae | Leptolobium elegans | 2 | 32 | 34 | 1 | 17 | 18 | 3 | 32 | 35 | 3 | 0 | 3 | 1 | 0 | 1 | 10 | 81 | 91 |
| Fabaceae | Leucochloron incuriale | 0 | 0 | 0 | 0 | 5 | 5 | 0 | 0 | 0 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 7 | 7 |
| Fabaceae | Machaerium acutifolium | 1 | 0 | 1 | 1 | 3 | 4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 3 | 5 |
| Fabaceae | Machaerium brasiliense | 1 | 0 | 1 | 1 | 5 | 6 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 3 | 5 | 8 |
| Fabaceae | Machaerium villosum | 0 | 0 | 0 | 0 | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3 |
| Fabaceae | Platypodium elegans | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 2 | 0 | 2 |
| Fabaceae | Pterogyne nitens | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 |
| Fabaceae | Senegalia polyphylla | 48 | 0 | 48 | 94 | 1 | 95 | 87 | 2 | 89 | 101 | 0 | 101 | 0 | 0 | 0 | 330 | 3 | 333 |
| Fabaceae | Senna pendula | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 1 |
| Fabaceae | Senna rugosa | 6 | 0 | 6 | 7 | 0 | 7 | 9 | 0 | 9 | 7 | 35 | 42 | 0 | 28 | 28 | 29 | 63 | 92 |
| Lamiaceae | Aegiphila verticillata | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Lamiaceae | Vitex polygama | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 2 | 14 | 0 | 14 | 14 | 2 | 16 |
| Lauraceae | Aiouea sellowiana | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 4 | 4 |
| Lauraceae | Endlicheria paniculata | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 2 | 4 | 0 | 0 | 0 | 2 | 2 | 4 |
| Lauraceae | Ocotea corymbosa | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 1 | 3 |
| Lauraceae | Ocotea puberula | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 |
| Lauraceae | Ocotea pulchella | 9 | 5 | 14 | 2 | 0 | 2 | 2 | 1 | 3 | 0 | 6 | 6 | 0 | 0 | 0 | 13 | 12 | 25 |
| Lauraceae | Ocotea velutina | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 1 |
| Lythraceae | Lafoensia pacari | 0 | 0 | 0 | 1 | 1 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 2 |
| Malpighiaceae | Byrsonima crassifolia | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Malpighiaceae | Byrsonima intermedia | 0 | 0 | 0 | 6 | 1 | 7 | 3 | 0 | 3 | 0 | 6 | 6 | 0 | 0 | 0 | 9 | 7 | 16 |
| Malvaceae | Eriotheca gracilipes | 1 | 1 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 2 |
| Melastomataceae | Miconia albicans | 13 | 12 | 25 | 2 | 6 | 8 | 0 | 3 | 3 | 9 | 4 | 13 | 7 | 0 | 7 | 31 | 25 | 56 |
| Melastomataceae | Miconia auricoma | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 11 | 11 | 1 | 3 | 4 | 2 | 14 | 16 |
| Melastomataceae | Miconia langsdorffii | 0 | 5 | 5 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 | 2 | 5 | 7 |
| Melastomataceae | Miconia ligustroides | 5 | 0 | 5 | 5 | 0 | 5 | 0 | 0 | 0 | 6 | 7 | 13 | 0 | 2 | 2 | 16 | 9 | 25 |
| Melastomataceae | Miconia minutiflora | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3 | 2 | 0 | 2 | 2 | 3 | 5 |
| Myrtaceae | Campomanesia pubescens | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 2 | 0 | 0 | 0 | 1 | 1 | 2 |
| Myrtaceae | Eugenia aurata | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 6 | 7 | 0 | 0 | 0 | 2 | 6 | 8 |
| Myrtaceae | Eugenia bimarginata | 1 | 0 | 1 | 1 | 2 | 3 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 3 | 2 | 5 |
| Myrtaceae | Eugenia hyemalis | 1 | 0 | 1 | 9 | 0 | 9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 10 | 0 | 10 |
| Myrtaceae | Eugenia myrcianthes | 0 | 0 | 0 | 2 | 1 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 1 | 3 |
| Myrtaceae | Eugenia punicifolia | 5 | 0 | 5 | 6 | 0 | 6 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 12 | 0 | 12 |
| Myrtaceae | Eugenia pyriformis | 1 | 0 | 1 | 6 | 0 | 6 | 0 | 0 | 0 | 6 | 0 | 6 | 0 | 0 | 0 | 13 | 0 | 13 |
| Myrtaceae | Eugenia uniflora | 0 | 0 | 0 | 5 | 0 | 5 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 5 | 0 | 5 |
| Myrtaceae | Myrcia | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 |
| Myrtaceae | Myrcia bella | 2 | 14 | 16 | 3 | 11 | 14 | 5 | 8 | 13 | 2 | 3 | 5 | 0 | 0 | 0 | 12 | 36 | 48 |
| Myrtaceae | Myrcia guianensis | 5 | 5 | 10 | 1 | 0 | 1 | 4 | 2 | 6 | 0 | 0 | 0 | 0 | 0 | 0 | 10 | 7 | 17 |
| Myrtaceae | Myrcia rufipes | 1 | 4 | 5 | 2 | 0 | 2 | 1 | 0 | 1 | 3 | 0 | 3 | 0 | 0 | 0 | 7 | 4 | 11 |
| Myrtaceae | Myrcia splendens | 4 | 0 | 4 | 4 | 0 | 4 | 0 | 1 | 1 | 4 | 3 | 7 | 0 | 0 | 0 | 12 | 4 | 16 |
| Myrtaceae | Myrcia tomentosa | 0 | 0 | 0 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 3 | 3 |
| Myrtaceae | Myrciaria floribunda | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 6 | 6 | 0 | 0 | 0 | 1 | 6 | 7 |
| Myrtaceae | Myrciaria glazioviana | 11 | 0 | 11 | 4 | 1 | 5 | 5 | 0 | 5 | 7 | 0 | 7 | 0 | 0 | 0 | 27 | 1 | 28 |
| Myrtaceae | Psidium grandifolium | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 2 | 0 | 2 |
| Myrtaceae | Psidium guineense | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 |
| Myrtaceae | Psidium oligospermum | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 |
| Nyctaginaceae | Guapira noxia | 0 | 4 | 4 | 0 | 1 | 1 | 0 | 0 | 0 | 3 | 0 | 3 | 0 | 0 | 0 | 3 | 5 | 8 |
| Ochnaceae | Ouratea spectabilis | 9 | 12 | 21 | 4 | 9 | 13 | 1 | 0 | 1 | 6 | 0 | 6 | 0 | 0 | 0 | 20 | 21 | 41 |
| Peraceae | Pera glabrata | 0 | 0 | 0 | 1 | 0 | 1 | 4 | 8 | 12 | 1 | 1 | 2 | 0 | 2 | 2 | 6 | 11 | 17 |
| Piperaceae | Piper mollicomum | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Polygalaceae | Bredemeyera floribunda | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Polygonaceae | Coccoloba cordata | 0 | 0 | 0 | 3 | 5 | 8 | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 5 | 5 | 10 |
| Polygonaceae | Triplaris americana | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Primulaceae | Myrsine coriacea | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 6 | 8 | 0 | 0 | 0 | 2 | 6 | 8 |
| Primulaceae | Myrsine gardneriana | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Primulaceae | Myrsine guianensis | 0 | 0 | 0 | 4 | 0 | 4 | 0 | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 4 | 3 | 7 |
| Primulaceae | Myrsine umbellata | 2 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 0 | 3 | 0 | 0 | 0 | 5 | 0 | 5 |
| Proteaceae | Roupala montana | 1 | 2 | 3 | 0 | 3 | 3 | 0 | 0 | 0 | 3 | 0 | 3 | 0 | 0 | 0 | 4 | 5 | 9 |
| Rubiaceae | Amaioua intermedia | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 4 | 0 | 4 | 18 | 0 | 18 | 23 | 0 | 23 |
| Rubiaceae | Palicourea sessilis | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 4 | 1 | 5 | 0 | 0 | 0 | 5 | 1 | 6 |
| Salicaceae | Casearia sylvestris | 36 | 0 | 36 | 76 | 0 | 76 | 60 | 0 | 60 | 15 | 0 | 15 | 2 | 0 | 2 | 189 | 0 | 189 |
| Sapindaceae | Matayba elaeagnoides | 0 | 0 | 0 | 3 | 0 | 3 | 2 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 5 | 0 | 5 |
| Sapindaceae | Sapindus saponaria | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 3 | 3 | 11 | 0 | 11 | 12 | 3 | 15 |
| Sapotaceae | Pouteria ramiflora | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 1 | 1 | 26 | 0 | 26 | 26 | 2 | 28 |
| Solanaceae | Cestrum corymbosum | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3 | 3 | 0 | 3 | 3 | 3 | 6 |
| Solanaceae | Solanum didymum | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 |
| Solanaceae | Solanum granuloso-leprosum | 44 | 3 | 47 | 79 | 6 | 85 | 77 | 1 | 78 | 55 | 3 | 58 | 2 | 7 | 9 | 257 | 20 | 277 |
| Solanaceae | Solanum sisymbriifolium | 166 | 0 | 166 | 72 | 0 | 72 | 80 | 0 | 80 | 98 | 0 | 98 | 0 | 0 | 0 | 416 | 0 | 416 |
| Urticaceae | Cecropia pachystachya | 1 | 0 | 1 | 2 | 0 | 2 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 9 | 9 | 4 | 9 | 13 |
| Verbenaceae | Lantana camara | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 0 | 0 | 0 | 12 | 12 | 1 | 0 | 1 | 3 | 12 | 15 |
| Vochysiaceae | Qualea grandiflora | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 0 | 0 | 222 | 0 | 222 | 224 | 0 | 224 |
| NA | Unclassified | 9 | 8 | 17 | 5 | 4 | 9 | 8 | 1 | 9 | 12 | 2 | 14 | 0 | 0 | 0 | 34 | 15 | 49 |
! Altogether
sampled_species <-readxl::read_xlsx("itatinga_sampled_species.xlsx")
sampled_species |>
dplyr::select(-`Seed dispersal`) |>
knitr::kable(format = "html")| family | scientificName | scientificNameAuthorship | Large_recruit | Small_recruit | Origin | Planting group | Successional group |
|---|---|---|---|---|---|---|---|
| Anacardiaceae | Astronium graveolens | Jacq. | TRUE | TRUE | planted | diversity | non-pioneer |
| Anacardiaceae | Schinus terebinthifolia | Raddi | FALSE | TRUE | recruit | NA | pioneer |
| Anacardiaceae | Myracrodruon urundeuva | Allem. | TRUE | FALSE | recruit | NA | non-pioneer |
| Anacardiaceae | Tapirira guianensis | Aubl. | TRUE | FALSE | recruit | NA | non-pioneer |
| Annonaceae | Annona cacans | Warm. | TRUE | TRUE | planted | diversity | non-pioneer |
| Annonaceae | Annona coriacea | Mart. | TRUE | TRUE | recruit | NA | non-pioneer |
| Annonaceae | Duguetia furfuracea | (A.St.-Hil.) Saff. | TRUE | TRUE | recruit | NA | non-pioneer |
| Apocynaceae | Aspidosperma tomentosum | Mart. | TRUE | TRUE | recruit | NA | non-pioneer |
| Araliaceae | Didymopanax vinosus | (Cham. & Schltdl.) Marchal | TRUE | TRUE | recruit | NA | non-pioneer |
| Arecaceae | Syagrus romanzoffiana | (Cham.) Glassman | TRUE | TRUE | recruit | NA | non-pioneer |
| Asteraceae | Piptocarpha rotundifolia | Baker | FALSE | TRUE | recruit | NA | pioneer |
| Asteraceae | Baccharis dracunculifolia | DC. | TRUE | TRUE | recruit | NA | pioneer |
| Asteraceae | Moquiniastrum barrosoae | (Cabrera) G.Sancho | TRUE | TRUE | recruit | NA | pioneer |
| Asteraceae | Moquiniastrum polymorphum | (Less.) G.Sancho | TRUE | TRUE | recruit | NA | pioneer |
| Asteraceae | Vernonanthura polyanthes | (Spreng.) A.J.Vega & Dematt. | TRUE | TRUE | recruit | NA | pioneer |
| Bignoniaceae | Handroanthus ochraceus | (Cham.) Mattos | TRUE | TRUE | planted | diversity | non-pioneer |
| Bignoniaceae | Handroanthus impetiginosus | (Mart. ex DC.) Mattos | FALSE | FALSE | planted | diversity | non-pioneer |
| Bignoniaceae | Jacaranda cuspidifolia | Mart. | FALSE | FALSE | planted | diversity | non-pioneer |
| Bignoniaceae | Cybistax antisyphilitica | (Mart.) Mart. | TRUE | FALSE | recruit | NA | non-pioneer |
| Bignoniaceae | Handroanthus albus | (Cham.) Mattos | TRUE | TRUE | recruit | NA | non-pioneer |
| Boraginaceae | Cordia ecalyculata | Vell. | FALSE | FALSE | planted | diversity | pioneer |
| Boraginaceae | Cordia sellowiana | Cham. | TRUE | FALSE | recruit | NA | pioneer |
| Burseraceae | Protium spruceanum | Engl. | TRUE | FALSE | recruit | NA | non-pioneer |
| Cannabaceae | Trema micrantha | (L.) Blume | TRUE | FALSE | recruit | NA | pioneer |
| Chrysobalanaceae | Couepia grandiflora | (Mart. & Zucc.) Benth. ex Hook.f. | TRUE | TRUE | recruit | NA | non-pioneer |
| Combretaceae | Terminalia argentea | Mart. | TRUE | FALSE | recruit | NA | non-pioneer |
| Euphorbiaceae | Croton floribundus | Spreng. | TRUE | TRUE | planted | filling | pioneer |
| Euphorbiaceae | Croton urucurana | Baill. | TRUE | TRUE | planted | filling | pioneer |
| Euphorbiaceae | Mabea fistulifera | Mart. | TRUE | TRUE | planted | filling | pioneer |
| Euphorbiaceae | Alchornea triplinervia | (Spreng.) Müll.Arg. | FALSE | TRUE | recruit | NA | non-pioneer |
| Euphorbiaceae | Sapium glandulosum | (L.) Morong | TRUE | TRUE | recruit | NA | non-pioneer |
| Fabaceae | Enterolobium contortisiliquum | (Vell.) Morong | FALSE | TRUE | planted | diversity | non-pioneer |
| Fabaceae | Erythrina mulungu | Mart. ex Benth. | FALSE | TRUE | planted | diversity | non-pioneer |
| Fabaceae | Platypodium elegans | Vogel | FALSE | TRUE | planted | diversity | non-pioneer |
| Fabaceae | Copaifera langsdorffii | Desf. | TRUE | TRUE | planted | diversity | non-pioneer |
| Fabaceae | Machaerium acutifolium | Vogel | TRUE | TRUE | planted | diversity | non-pioneer |
| Fabaceae | Pterogyne nitens | Tul. | TRUE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Anadenanthera colubrina var. cebil | (Griseb.) Altschul | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Anadenanthera peregrina | (L.) Speg. | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Apuleia leiocarpa | (Vogel) J.F.Macbr. | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Cassia fistula | L. | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Hymenaea courbaril | L. | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Myroxylon peruiferum | L.f. | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Parapiptadenia rigida | (Benth.) Brenan | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Peltophorum dubium | (Spreng.) Taub. | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Senna macranthera | (DC. ex Collad.) H.S.Irwin & Barneby | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Inga vera | Willd. | FALSE | TRUE | planted | filling | non-pioneer |
| Fabaceae | Senegalia polyphylla | (DC.) Britton & Rose | TRUE | TRUE | planted | filling | pioneer |
| Fabaceae | Inga edulis | Mart. | FALSE | FALSE | planted | filling | non-pioneer |
| Fabaceae | Andira fraxinifolia | Benth. | FALSE | TRUE | recruit | NA | non-pioneer |
| Fabaceae | Bauhinia rufa | (Bong.) Steud. | FALSE | TRUE | recruit | NA | non-pioneer |
| Fabaceae | Chamaecrista comosa | E.Mey. | TRUE | FALSE | recruit | NA | non-pioneer |
| Fabaceae | Dalbergia frutescens | (Vell.) Britton | TRUE | TRUE | recruit | NA | non-pioneer |
| Fabaceae | Dalbergia miscolobium | Benth. | TRUE | TRUE | recruit | NA | non-pioneer |
| Fabaceae | Hymenaea stigonocarpa | Mart. ex Hayne | TRUE | FALSE | recruit | NA | non-pioneer |
| Fabaceae | Leptolobium dasycarpum | Vogel | TRUE | FALSE | recruit | NA | non-pioneer |
| Fabaceae | Leptolobium elegans | Vogel | TRUE | TRUE | recruit | NA | non-pioneer |
| Fabaceae | Leucochloron incuriale | (Vell.) Barneby & J.W.Grimes | TRUE | FALSE | recruit | NA | non-pioneer |
| Fabaceae | Machaerium brasiliense | Vogel | TRUE | TRUE | recruit | NA | non-pioneer |
| Fabaceae | Machaerium villosum | Vogel | TRUE | FALSE | recruit | NA | non-pioneer |
| Fabaceae | Senna pendula | (Humb. & Bonpl. ex Willd.) H.S.Irwin & Barneby | TRUE | FALSE | recruit | NA | pioneer |
| Fabaceae | Senna rugosa | (G.Don) H.S.Irwin & Barneby | TRUE | TRUE | recruit | NA | pioneer |
| Lamiaceae | Aegiphila verticillata | Vell. | FALSE | TRUE | recruit | NA | pioneer |
| Lamiaceae | Vitex polygama | Cham. | TRUE | TRUE | recruit | NA | non-pioneer |
| Lauraceae | Ocotea puberula | (Rich.) Nees | FALSE | TRUE | recruit | NA | non-pioneer |
| Lauraceae | Aiouea sellowiana | (Nees & Mart.) R.Rohde | TRUE | FALSE | recruit | NA | non-pioneer |
| Lauraceae | Endlicheria paniculata | (Spreng.) J.F.Macbr. | TRUE | TRUE | recruit | NA | non-pioneer |
| Lauraceae | Ocotea corymbosa | (Meisn.) Mez | TRUE | TRUE | recruit | NA | non-pioneer |
| Lauraceae | Ocotea pulchella | (Nees & Mart.) Mez | TRUE | TRUE | recruit | NA | non-pioneer |
| Lauraceae | Ocotea velutina | (Nees) Mart. ex B.D.Jacks. | TRUE | FALSE | recruit | NA | non-pioneer |
| Lecythidaceae | Cariniana estrellensis | (Raddi) Kuntze | FALSE | FALSE | planted | diversity | non-pioneer |
| Lythraceae | Lafoensia pacari | A.St.-Hil. | TRUE | TRUE | planted | diversity | non-pioneer |
| Malpighiaceae | Byrsonima crassifolia | Kunth | FALSE | TRUE | recruit | NA | non-pioneer |
| Malpighiaceae | Byrsonima intermedia | A.Juss. | TRUE | TRUE | recruit | NA | non-pioneer |
| Malvaceae | Ceiba speciosa | (A.St.-Hil., A.Juss. & Cambess.) Ravenna | FALSE | FALSE | planted | diversity | non-pioneer |
| Malvaceae | Heliocarpus popayanensis | Kunth | FALSE | FALSE | planted | diversity | pioneer |
| Malvaceae | Luehea divaricata | Mart. | FALSE | FALSE | planted | diversity | non-pioneer |
| Malvaceae | Guazuma ulmifolia | Lam. | FALSE | FALSE | planted | filling | pioneer |
| Malvaceae | Eriotheca gracilipes | (K.Schum.) A.Robyns | TRUE | TRUE | recruit | NA | non-pioneer |
| Melastomataceae | Miconia albicans | (Sw.) Steud. | TRUE | TRUE | recruit | NA | pioneer |
| Melastomataceae | Miconia auricoma | (Spring ex Mart.) R.Goldenb. | TRUE | TRUE | recruit | NA | pioneer |
| Melastomataceae | Miconia langsdorffii | Cogn. | TRUE | TRUE | recruit | NA | pioneer |
| Melastomataceae | Miconia ligustroides | Naudin | TRUE | TRUE | recruit | NA | pioneer |
| Melastomataceae | Miconia minutiflora | (Bonpl.) DC. | TRUE | TRUE | recruit | NA | pioneer |
| Meliaceae | Cedrela fissilis | Vell. | FALSE | FALSE | planted | diversity | non-pioneer |
| Moraceae | Ficus guaranitica | Chodat | FALSE | FALSE | planted | diversity | non-pioneer |
| Myrtaceae | Eugenia uniflora | L. | FALSE | TRUE | planted | diversity | non-pioneer |
| Myrtaceae | Psidium myrtoides | O.Berg | FALSE | FALSE | planted | diversity | non-pioneer |
| Myrtaceae | Eugenia hyemalis | Cambess. | FALSE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Eugenia punicifolia | (Kunth) DC. | FALSE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Eugenia pyriformis | Cambess. | FALSE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Myrcia | DC. | FALSE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Psidium grandifolium | Mart. ex DC. | FALSE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Psidium guineense | Sw. | FALSE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Campomanesia pubescens | O.Berg | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Eugenia aurata | O.Berg | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Eugenia bimarginata | DC. | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Eugenia myrcianthes | Nied. | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Myrcia bella | Cambess. | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Myrcia guianensis | DC. | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Myrcia rufipes | DC. | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Myrcia splendens | DC. | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Myrcia tomentosa | DC. | TRUE | FALSE | recruit | NA | non-pioneer |
| Myrtaceae | Myrciaria floribunda | O.Berg | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Myrciaria glazioviana | (Kiaersk.) G.M.Barroso ex Sobral | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Psidium oligospermum | Mart. ex DC. | TRUE | FALSE | recruit | NA | non-pioneer |
| Nyctaginaceae | Guapira noxia | (Netto) Lundell | TRUE | TRUE | recruit | NA | non-pioneer |
| Ochnaceae | Ouratea spectabilis | Engl. | TRUE | TRUE | recruit | NA | non-pioneer |
| Peraceae | Pera glabrata | (Schott) Poepp. ex Baill. | TRUE | TRUE | recruit | NA | non-pioneer |
| Petiveriaceae | Gallesia integrifolia | (Spreng.) Harms | FALSE | FALSE | planted | diversity | non-pioneer |
| Piperaceae | Piper mollicomum | Kunth | FALSE | TRUE | recruit | NA | non-pioneer |
| Polygalaceae | Bredemeyera floribunda | Willd. | FALSE | TRUE | recruit | NA | non-pioneer |
| Polygonaceae | Triplaris americana | L. | FALSE | TRUE | recruit | NA | non-pioneer |
| Polygonaceae | Coccoloba cordata | Cham. | TRUE | TRUE | recruit | NA | non-pioneer |
| Primulaceae | Myrsine gardneriana | A.DC. | FALSE | TRUE | recruit | NA | non-pioneer |
| Primulaceae | Myrsine gardneriana | A.DC. | FALSE | TRUE | recruit | NA | non-pioneer |
| Primulaceae | Myrsine umbellata | Mart. | FALSE | TRUE | recruit | NA | non-pioneer |
| Primulaceae | Myrsine coriacea | (Sw.) R.Br. ex Roem. & Schult. | TRUE | TRUE | recruit | NA | non-pioneer |
| Primulaceae | Myrsine guianensis | (Aubl.) Kuntze | TRUE | TRUE | recruit | NA | non-pioneer |
| Proteaceae | Roupala montana | Aubl. | TRUE | TRUE | recruit | NA | non-pioneer |
| Rosaceae | Prunus myrtifolia | (L.) Urb. | FALSE | FALSE | planted | diversity | non-pioneer |
| Rubiaceae | Genipa americana | L. | FALSE | FALSE | planted | diversity | non-pioneer |
| Rubiaceae | Amaioua intermedia | Mart. | FALSE | TRUE | recruit | NA | non-pioneer |
| Rubiaceae | Palicourea sessilis | (Vell.) C.M.Taylor | TRUE | TRUE | recruit | NA | non-pioneer |
| Salicaceae | Casearia sylvestris | Sw. | FALSE | TRUE | recruit | NA | pioneer |
| Sapindaceae | Matayba elaeagnoides | Radlk. | FALSE | TRUE | recruit | NA | non-pioneer |
| Sapindaceae | Sapindus saponaria | L. | TRUE | TRUE | recruit | NA | non-pioneer |
| Sapotaceae | Pouteria ramiflora | Radlk. | TRUE | TRUE | recruit | NA | non-pioneer |
| Solanaceae | Solanum didymum | Dunal | FALSE | TRUE | recruit | NA | pioneer |
| Solanaceae | Solanum sisymbriifolium | Lam. | FALSE | TRUE | recruit | NA | pioneer |
| Solanaceae | Cestrum corymbosum | Schltdl. | TRUE | FALSE | recruit | NA | non-pioneer |
| Solanaceae | Solanum granuloso-leprosum | Dunal | TRUE | TRUE | recruit | NA | pioneer |
| Unclassified | Unclassified | NA | TRUE | TRUE | recruit | NA | NA |
| Urticaceae | Cecropia pachystachya | Trécul | TRUE | TRUE | planted | diversity | pioneer |
| Verbenaceae | Citharexylum myrianthum | Cham. | FALSE | FALSE | planted | diversity | non-pioneer |
| Verbenaceae | Lantana camara | L. | TRUE | TRUE | recruit | NA | pioneer |
| Vochysiaceae | Qualea grandiflora | Mart. | FALSE | TRUE | recruit | NA | non-pioneer |
Estimation in Mg.ha-1
576 = 12m x 12m x 4 plots
AGB_data0 <- trees_clean |>
dplyr::group_by(ano, tratamento, bloco, plantio)|>
dplyr::summarise(AGB_sum = sum(AGB_ferez, na.rm=T)/1000, .groups = "drop") |> # writexl::write_xlsx("AGB_plot.xlsx")
dplyr::mutate(
AGB_Mgha = dplyr::case_when(
tratamento == "100" ~ (AGB_sum / 576) * 10000,
tratamento %in% c("F50", "N50") & plantio == "CP" ~ (AGB_sum / 576) * 5000,
tratamento %in% c("F50", "N50") & plantio == "SP" ~ (AGB_sum / 576) * 5000,
tratamento %in% c("F25", "N25") & plantio == "CP" ~ (AGB_sum / 576) * 2500,
tratamento %in% c("F25", "N25") & plantio == "SP" ~ (AGB_sum / 576) * 7500),
Tratamento = ordered(tratamento, levels = ordem,
labels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%")))
AGB_data <- AGB_data0 |>
dplyr::group_by(ano, tratamento, bloco)|>
dplyr::summarise(AGB_Mgha = sum(AGB_Mgha), .groups = "drop") |>
dplyr::mutate(
Ano = ano-2021,
ANO = factor(ano),
scaled_AGB = scale(AGB_Mgha))Comparing models
AGB_model0 <- lme4::lmer(scaled_AGB ~ 0 + (1|bloco), data = AGB_data)
AGB_model1 <- lme4::lmer(scaled_AGB ~ ANO + (1|bloco), data = AGB_data)
AGB_model2 <- lme4::lmer(scaled_AGB ~ ANO + bloco + (1|bloco), data = AGB_data) Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Hessian is numerically singular: parameters are not uniquely determined
AGB_model3 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + bloco + (1|bloco), data = AGB_data) # consider rescaling
AGB_model4 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + bloco + (1|bloco), data = AGB_data) # isSingular
AGB_model5 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + (1|bloco), data = AGB_data) # best model
AGB_model6 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1|bloco), data = AGB_data)
AGB_model7 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1+Ano|bloco), data = AGB_data) # isSingularboundary (singular) fit: see help('isSingular')
AIC(AGB_model0, AGB_model1, AGB_model2, AGB_model3, AGB_model4, AGB_model5,AGB_model6, AGB_model7 ) df AIC
AGB_model0 2 144.60044
AGB_model1 4 141.90535
AGB_model2 8 144.14482
AGB_model3 16 93.47468
AGB_model4 12 88.21817
AGB_model5 8 85.97871
AGB_model6 12 91.23521
AGB_model7 14 90.49850
Exploring residuals
par(mfrow=c(1,2))
DHARMa::plotResiduals(AGB_model5)
DHARMa::plotQQunif(AGB_model5)One outlier
Marginal \(R^2\) / Conditional \(R^2\)
MuMIn::r.squaredGLMM(AGB_model5) # without interaction R2m R2c
[1,] 0.6955092 0.8374889
0.696/ 0.837[1] 0.8315412
MuMIn::r.squaredGLMM(AGB_model6) # with interaction R2m R2c
[1,] 0.701532 0.8419478
0.702/ 0.842[1] 0.8337292
0.702/ 0.842 - variance explained by fixed effects and by the entire model (around r0.702/ 0.842 % is explained by fixed effects)
car::Anova(AGB_model6)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: scaled_AGB
Chisq Df Pr(>Chisq)
tratamento 165.4367 4 < 2.2e-16 ***
ANO 47.4491 1 5.645e-12 ***
tratamento:ANO 4.6061 4 0.3302
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction is no significant at all (\(p-value \sim 0.33\))
car::Anova(AGB_model5)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: scaled_AGB
Chisq Df Pr(>Chisq)
tratamento 162.968 4 < 2.2e-16 ***
ANO 46.741 1 8.102e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Summary
sjPlot::tab_model(AGB_model5)| scaled AGB | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 0.96 | 0.51 – 1.42 | <0.001 |
| tratamento [F25] | -1.83 | -2.20 – -1.46 | <0.001 |
| tratamento [F50] | -1.35 | -1.72 – -0.98 | <0.001 |
| tratamento [N25] | -2.19 | -2.56 – -1.82 | <0.001 |
| tratamento [N50] | -1.44 | -1.82 – -1.07 | <0.001 |
| ANO [2024] | 0.80 | 0.56 – 1.03 | <0.001 |
| Random Effects | |||
| σ2 | 0.17 | ||
| τ00 bloco | 0.15 | ||
| ICC | 0.47 | ||
| N bloco | 5 | ||
| Observations | 50 | ||
| Marginal R2 / Conditional R2 | 0.696 / 0.837 | ||
broom.mixed::tidy(AGB_model5) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 8 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 0.965 0.224 4.31
2 fixed <NA> tratamentoF25 -1.83 0.184 -9.92
3 fixed <NA> tratamentoF50 -1.35 0.184 -7.34
4 fixed <NA> tratamentoN25 -2.19 0.184 -11.9
5 fixed <NA> tratamentoN50 -1.44 0.184 -7.83
6 fixed <NA> ANO2024 0.797 0.117 6.84
7 ran_pars bloco sd__(Intercept) 0.385 NA NA
8 ran_pars Residual sd__Observation 0.412 NA NA
Assessing comparisons:
emmeans::emmeans(AGB_model6, pairwise ~ tratamento | ANO) |> multcomp::cld()ANO = 2022:
tratamento emmean SE df lower.CL upper.CL .group
N25 -1.028 0.251 13.4 -1.5696 -0.4868 1
F25 -0.804 0.251 13.4 -1.3450 -0.2622 1
N50 -0.527 0.251 13.4 -1.0680 0.0148 1
F50 -0.423 0.251 13.4 -0.9645 0.1183 1
100 0.790 0.251 13.4 0.2484 1.3312 2
ANO = 2024:
tratamento emmean SE df lower.CL upper.CL .group
N25 -0.629 0.251 13.4 -1.1701 -0.0873 1
F25 -0.125 0.251 13.4 -0.6663 0.4165 12
N50 0.366 0.251 13.4 -0.1756 0.9072 2
F50 0.443 0.251 13.4 -0.0981 0.9847 2
100 1.936 0.251 13.4 1.3949 2.4777 3
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
AGB_model6a <- lme4::lmer(scaled_AGB ~ tratamento*Ano + (1|bloco), data = AGB_data)
AIC(AGB_model6a, AGB_model6) df AIC
AGB_model6a 12 98.16669
AGB_model6 12 91.23521
emmeans::emtrends(AGB_model6a, pairwise ~ tratamento, var = "Ano", poly.trend = TRUE) |> multcomp::cld() tratamento Ano.trend SE df lower.CL upper.CL .group
N25 0.200 0.129 36 -0.0625 0.462 1
F25 0.339 0.129 36 0.0771 0.602 1
F50 0.433 0.129 36 0.1709 0.695 1
N50 0.446 0.129 36 0.1840 0.708 1
100 0.573 0.129 36 0.3110 0.836 1
Results are averaged over the levels of: Ano
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Preparing data:
Comparison only inside plots over the years
ordem <-c("N25", "F25", "N50", "F50", "100")
AGB_inside0 <- trees_clean |>
dplyr::group_by(ano, tratamento, bloco, parcela, plantio)|>
dplyr::summarise(
AGB_sum = sum(AGB_ferez, na.rm=T),
.groups = "drop") |>
dplyr::mutate(
tratamento = ordered(tratamento, levels = ordem),
Parcela = paste(tratamento, bloco, parcela, sep="_"),
Ano = ano-2021) |>
dplyr::filter(plantio == "CP")View
AGB_inside0 |>
ggplot(aes(x = factor(ano), y = log(AGB_sum), group = Parcela, colour = parcela)) +
geom_point() +
geom_line()+facet_grid(bloco~tratamento)+
theme_minimal()Plots with decreasing :
deltas <- AGB_inside0 |>
dplyr::select(ano, Parcela, AGB_sum)|>
tidyr::pivot_wider(names_from = ano, names_prefix = "ano_", values_from = AGB_sum)|>
dplyr::mutate(delta_AGB = (ano_2024/ano_2022))|>
dplyr::arrange(delta_AGB)
deltas |> head(10)# A tibble: 10 × 4
Parcela ano_2022 ano_2024 delta_AGB
<chr> <dbl> <dbl> <dbl>
1 100_4_4 422. 216. 0.511
2 F50_5_3 184. 128. 0.698
3 N50_5_4 324. 228. 0.701
4 F25_5_1 269. 198. 0.736
5 F50_1_1 152. 141. 0.930
6 N25_3_4 165. 157. 0.954
7 100_5_3 230. 220. 0.957
8 F50_3_4 208. 208. 0.998
9 N50_1_4 485. 493. 1.02
10 F50_2_4 549. 577. 1.05
deltas |> tail(10)# A tibble: 10 × 4
Parcela ano_2022 ano_2024 delta_AGB
<chr> <dbl> <dbl> <dbl>
1 F50_2_1 216. 758. 3.51
2 N50_5_1 170. 611. 3.58
3 N50_5_3 153. 554. 3.62
4 N25_3_2 223. 846. 3.79
5 F50_2_3 77.5 313. 4.03
6 F50_2_2 82.8 367. 4.43
7 F50_3_3 175. 799. 4.56
8 F25_3_3 114. 557. 4.88
9 N50_4_2 136. 784. 5.75
10 N50_2_4 88.6 690. 7.79
How much % is planted or recruit in AGB?
trees_clean |>
dplyr::filter(plantio == "CP") |>
#dplyr::filter(ano == 2024) |>
dplyr::group_by(ano, tratamento, bloco)|>
dplyr::summarise(
AGB_sum = sum(AGB_ferez, na.rm=T)/1000,
AGB_planted = sum(AGB_ferez[RP == "planted"], na.rm = TRUE)/1000,
prop_AGB_planted = (AGB_planted / AGB_sum)*100,
.groups = "drop") |> #dplyr::group_by(tratamento,plantio)|>
dplyr::summarise(mean(prop_AGB_planted),
sd(prop_AGB_planted)/sqrt(dplyr::n()))# A tibble: 1 × 2
`mean(prop_AGB_planted)` `sd(prop_AGB_planted)/sqrt(dplyr::n())`
<dbl> <dbl>
1 96.9 0.758
Should I treat them with a random factor (1|bloco)?
AGB_inside <- trees_clean |>
dplyr::filter(plantio == "CP") |>
#dplyr::filter(ano == 2024) |>
dplyr::group_by(ano, tratamento, bloco)|>
dplyr::summarise(AGB_sum = sum(AGB_ferez, na.rm=T)/1000, .groups = "drop") |>
dplyr::mutate(
Tratamento = ordered(tratamento, levels = ordem),
Ano = ano-2021,
ANO = ordered(ano),
scaled_AGB = scale(AGB_sum)) Comparing models
AGB_inside_m0 <- lme4::lmer(scaled_AGB ~ 0 + (1|bloco), data = AGB_inside)
AGB_inside_m1 <- lme4::lmer(scaled_AGB ~ ANO + (1|bloco), data = AGB_inside)
AGB_inside_m2 <- lme4::lmer(scaled_AGB ~ ANO + bloco + (1|bloco), data = AGB_inside) # best model Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Hessian is numerically singular: parameters are not uniquely determined
AGB_inside_m3 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + bloco + (1|bloco), data = AGB_inside)
AGB_inside_m4 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + bloco + (1|bloco), data = AGB_inside)
AGB_inside_m5 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + (1|bloco), data = AGB_inside)
AGB_inside_m6 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1|bloco), data = AGB_inside)
AGB_inside_m7 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1+Ano|bloco), data = AGB_inside) # isSingularboundary (singular) fit: see help('isSingular')
AIC(AGB_inside_m0, AGB_inside_m1, AGB_inside_m2, AGB_inside_m3, AGB_inside_m4,AGB_inside_m5, AGB_inside_m6, AGB_inside_m7) df AIC
AGB_inside_m0 2 138.7732
AGB_inside_m1 4 117.1008
AGB_inside_m2 8 116.0631
AGB_inside_m3 16 122.2545
AGB_inside_m4 12 115.0812
AGB_inside_m5 8 116.1189
AGB_inside_m6 12 123.2922
AGB_inside_m7 14 123.2583
Exploring residuals
par(mfrow=c(1,2))
DHARMa::plotResiduals(AGB_inside_m6)
DHARMa::plotQQunif(AGB_inside_m6)Marginal \(R^2\) / Conditional \(R^2\)
MuMIn::r.squaredGLMM(AGB_inside_m6) R2m R2c
[1,] 0.396267 0.699652
0.396/ 0.6996[1] 0.5660377
0.396/ 0.6996 variance explained by fixed effects and by the entire model (around 0.5660377 % is explained by fixed effects)
car::Anova(AGB_inside_m6)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: scaled_AGB
Chisq Df Pr(>Chisq)
tratamento 14.3671 4 0.006211 **
ANO 47.6142 1 5.189e-12 ***
tratamento:ANO 2.6673 4 0.614946
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction is no significant at all (\(p-value \sim 0.85\))
sjPlot::tab_model(AGB_inside_m5)| scaled AGB | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 0.02 | -0.62 – 0.66 | 0.949 |
| tratamento [F25] | 0.50 | -0.02 – 1.01 | 0.059 |
| tratamento [F50] | 0.00 | -0.51 – 0.52 | 0.987 |
| tratamento [N25] | -0.47 | -0.99 – 0.04 | 0.070 |
| tratamento [N50] | -0.13 | -0.64 – 0.39 | 0.622 |
| ANO [linear] | 0.80 | 0.57 – 1.03 | <0.001 |
| Random Effects | |||
| σ2 | 0.33 | ||
| τ00 bloco | 0.34 | ||
| ICC | 0.51 | ||
| N bloco | 5 | ||
| Observations | 50 | ||
| Marginal R2 / Conditional R2 | 0.390 / 0.702 | ||
broom.mixed::tidy(AGB_inside_m5) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 8 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 0.02 0.318 0.064
2 fixed <NA> tratamentoF25 0.496 0.256 1.94
3 fixed <NA> tratamentoF50 0.004 0.256 0.016
4 fixed <NA> tratamentoN25 -0.475 0.256 -1.86
5 fixed <NA> tratamentoN50 -0.127 0.256 -0.497
6 fixed <NA> ANO.L 0.802 0.114 7.02
7 ran_pars bloco sd__(Intercept) 0.585 NA NA
8 ran_pars Residual sd__Observation 0.571 NA NA
Assessing comparisons:
emmeans::emmeans(AGB_inside_m6, pairwise ~ tratamento |ANO) |> multcomp::cld()ANO = 2022:
tratamento emmean SE df lower.CL upper.CL .group
N25 -0.8983 0.369 12.2 -1.700 -0.0969 1
N50 -0.7449 0.369 12.2 -1.546 0.0564 1
F50 -0.5376 0.369 12.2 -1.339 0.2637 1
100 -0.3908 0.369 12.2 -1.192 0.4105 1
F25 -0.2642 0.369 12.2 -1.066 0.5371 1
ANO = 2024:
tratamento emmean SE df lower.CL upper.CL .group
N25 -0.0107 0.369 12.2 -0.812 0.7906 1
100 0.4316 0.369 12.2 -0.370 1.2329 12
N50 0.5318 0.369 12.2 -0.269 1.3332 12
F50 0.5866 0.369 12.2 -0.215 1.3879 12
F25 1.2965 0.369 12.2 0.495 2.0979 2
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
AGB_inside_m6a <- lme4::lmer(scaled_AGB ~ tratamento*Ano + (1|bloco), data = AGB_inside)
AIC(AGB_inside_m6a, AGB_inside_m6) df AIC
AGB_inside_m6a 12 126.7579
AGB_inside_m6 12 123.2922
emmeans::emtrends(AGB_inside_m6a, pairwise ~ tratamento, var = "Ano", poly.trend = TRUE) |> multcomp::cld() tratamento Ano.trend SE df lower.CL upper.CL .group
100 0.411 0.184 36 0.0385 0.784 1
N25 0.444 0.184 36 0.0710 0.817 1
F50 0.562 0.184 36 0.1894 0.935 1
N50 0.638 0.184 36 0.2656 1.011 1
F25 0.780 0.184 36 0.4076 1.153 1
Results are averaged over the levels of: Ano
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Preparing data:
Comparison only outside plots over the years
AGB_outside0 <- trees_clean |>
dplyr::group_by(ano, tratamento, bloco, parcela, plantio)|>
dplyr::summarise(
AGB_sum = sum(AGB_ferez, na.rm=T),
.groups = "drop") |>
dplyr::mutate(
tratamento = ordered(tratamento, levels = ordem),
Parcela = paste(tratamento, bloco, parcela, sep="_"),
Ano = ano-2021) |>
dplyr::filter(plantio == "SP")View
AGB_outside0 |>
ggplot(aes(x = factor(ano), y = log(AGB_sum), group = Parcela, colour = parcela)) +
geom_point() +
geom_line()+facet_grid(bloco~tratamento)+
theme_minimal()it is a mess
Tratar os valores como (1|bloco)?
AGB_outside <- trees_clean |>
dplyr::filter(plantio == "SP") |>
#dplyr::filter(ano == 2024) |>
dplyr::group_by(ano, tratamento, bloco)|>
dplyr::summarise(AGB_sum = sum(AGB_ferez, na.rm=T)/1000, .groups = "drop") |>
dplyr::mutate(
Tratamento = ordered(tratamento, levels = ordem),
Ano = ano-2021,
ANO = ordered(ano),
scaled_AGB = scale(AGB_sum)) Comparing models
AGB_outside_m0 <- lme4::lmer(scaled_AGB ~ 0 + (1|bloco), data = AGB_outside) # best modelboundary (singular) fit: see help('isSingular')
AGB_outside_m1 <- lme4::lmer(scaled_AGB ~ ANO + (1|bloco), data = AGB_outside)
AGB_outside_m2 <- lme4::lmer(scaled_AGB ~ ANO + bloco + (1|bloco), data = AGB_outside) # best model
AGB_outside_m3 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + bloco + (1|bloco), data = AGB_outside)
AGB_outside_m4 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + bloco + (1|bloco), data = AGB_outside) Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradient
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Hessian is numerically singular: parameters are not uniquely determined
AGB_outside_m5 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + (1|bloco), data = AGB_outside)
AGB_outside_m6 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1|bloco), data = AGB_outside)
AGB_outside_m7 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1+Ano|bloco), data = AGB_outside) # isSingularboundary (singular) fit: see help('isSingular')
AIC(AGB_outside_m0, AGB_outside_m1, AGB_outside_m2, AGB_outside_m3,
AGB_outside_m4,AGB_outside_m5, AGB_outside_m6, AGB_outside_m7) df AIC
AGB_outside_m0 2 116.5024
AGB_outside_m1 4 117.0739
AGB_outside_m2 8 120.3459
AGB_outside_m3 14 128.0330
AGB_outside_m4 11 126.0844
AGB_outside_m5 7 122.8125
AGB_outside_m6 10 124.7610
AGB_outside_m7 12 126.6309
Exploring residuals
par(mfrow=c(1,2))
DHARMa::plotResiduals(AGB_outside_m1)
DHARMa::plotQQunif(AGB_outside_m0)Marginal \(R^2\) / Conditional \(R^2\)
MuMIn::r.squaredGLMM(AGB_outside_m1) R2m R2c
[1,] 0.1505529 0.1714753
0.1505529 /0.1714753[1] 0.8779859
0.1505529 /0.1714753 variance explained by fixed effects and by the entire model (around 0.8779859 % is explained by fixed effects)
car::Anova(AGB_outside_m5)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: scaled_AGB
Chisq Df Pr(>Chisq)
tratamento 0.6275 3 0.89012
ANO 6.5922 1 0.01024 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction is no significant at all (\(p-value \sim 0.85\))
sjPlot::tab_model(AGB_outside_m5)| scaled AGB | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 0.20 | -0.42 – 0.83 | 0.514 |
| tratamento [F50] | -0.30 | -1.17 – 0.57 | 0.493 |
| tratamento [N25] | -0.23 | -1.10 – 0.64 | 0.597 |
| tratamento [N50] | -0.29 | -1.15 – 0.58 | 0.509 |
| ANO [linear] | 0.55 | 0.11 – 0.98 | 0.015 |
| Random Effects | |||
| σ2 | 0.91 | ||
| τ00 bloco | 0.01 | ||
| ICC | 0.01 | ||
| N bloco | 5 | ||
| Observations | 40 | ||
| Marginal R2 / Conditional R2 | 0.154 / 0.167 | ||
broom.mixed::tidy(AGB_outside_m5) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 7 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 0.202 0.307 0.66
2 fixed <NA> tratamentoF50 -0.296 0.427 -0.694
3 fixed <NA> tratamentoN25 -0.228 0.427 -0.534
4 fixed <NA> tratamentoN50 -0.285 0.427 -0.667
5 fixed <NA> ANO.L 0.549 0.214 2.57
6 ran_pars bloco sd__(Intercept) 0.116 NA NA
7 ran_pars Residual sd__Observation 0.956 NA NA
Assessing comparisons:
emmeans::emmeans(AGB_outside_m6, pairwise ~ tratamento |ANO) |> multcomp::cld()ANO = 2022:
tratamento emmean SE df lower.CL upper.CL .group
F50 -0.6400 0.436 32 -1.5273 0.247 1
N25 -0.4221 0.436 32 -1.3094 0.465 1
F25 -0.3904 0.436 32 -1.2776 0.497 1
N50 -0.0995 0.436 32 -0.9868 0.788 1
ANO = 2024:
tratamento emmean SE df lower.CL upper.CL .group
N50 -0.0659 0.436 32 -0.9531 0.821 1
N25 0.3708 0.436 32 -0.5165 1.258 1
F50 0.4519 0.436 32 -0.4354 1.339 1
F25 0.7952 0.436 32 -0.0921 1.682 1
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 4 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Only inside
AGB_inside |>
dplyr::group_by(ano, tratamento) |>
dplyr::summarise(
Mean = round(mean(AGB_sum),2),
SE = round(sd(AGB_sum)/sqrt(dplyr::n()), 2),
.groups="drop") |>
dplyr::mutate(values = paste(Mean, "pm", SE)) |>
dplyr::select(-c(Mean, SE)) |>
tidyr::pivot_wider(names_from = tratamento, values_from = values)# A tibble: 2 × 6
ano `100` F25 F50 N25 N50
<dbl> <chr> <chr> <chr> <chr> <chr>
1 2022 1.23 pm 0.15 1.32 pm 0.2 1.12 pm 0.26 0.85 pm 0.2 0.96 pm 0.12
2 2024 1.85 pm 0.18 2.5 pm 0.37 1.96 pm 0.46 1.52 pm 0.37 1.92 pm 0.24
Increment per year:
AGB_inside |>
dplyr::group_by(tratamento) |>
dplyr::summarise(
variacao_media = mean((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2),
desvio_padrao = sd((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2),
erro_padrao = sd((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2) / sqrt(dplyr::n()),
min_variacao = min((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2),
max_variacao = max((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2),
.groups="drop")|>
dplyr::mutate(slope = paste(round(variacao_media,2), "pm", round(erro_padrao,2)))|>
dplyr::select(tratamento, slope)# A tibble: 5 × 2
tratamento slope
<chr> <chr>
1 100 0.31 pm 0.03
2 F25 0.59 pm 0.08
3 F50 0.42 pm 0.08
4 N25 0.33 pm 0.06
5 N50 0.48 pm 0.1
Total Aboveground biomass
tabS1_agb_total <- AGB_data |>
dplyr::group_by(ano, tratamento)|>
dplyr::summarise(Mean = round(mean(AGB_Mgha),2),
SE = round(sd(AGB_Mgha)/sqrt(dplyr::n()), 2),
.groups="drop")
tabS1_agb_total |>
dplyr::mutate(values = paste(Mean, "pm", SE)) |>
dplyr::select(-c(Mean, SE)) |>
tidyr::pivot_wider(names_from = tratamento, values_from = values)# A tibble: 2 × 6
ano `100` F25 F50 N25 N50
<dbl> <chr> <chr> <chr> <chr> <chr>
1 2022 21.35 pm 2.64 6.47 pm 0.92 10.03 pm 2.29 4.37 pm 1.07 9.06 pm 0.95
2 2024 32.06 pm 3.11 12.81 pm 2.12 18.12 pm 4.4 8.11 pm 1.64 17.39 pm 1.95
Increment per year:
AGB_data |>
dplyr::group_by(tratamento) %>%
dplyr::summarise(
variacao_media = mean((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2),
desvio_padrao = sd((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2),
erro_padrao = sd((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2) / sqrt(dplyr::n()),
min_variacao = min((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2),
max_variacao = max((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2),
.groups="drop")|>
dplyr::mutate(slope = paste(round(variacao_media,2), "pm", round(erro_padrao,2)))|>
dplyr::select(tratamento, slope)# A tibble: 5 × 2
tratamento slope
<chr> <chr>
1 100 5.35 pm 0.59
2 F25 3.17 pm 0.51
3 F50 4.05 pm 0.82
4 N25 1.87 pm 0.35
5 N50 4.17 pm 0.9
Percentuals
agb_plot_data <-AGB_data0 |>
dplyr::group_by(ano, Tratamento, plantio, bloco)|>
dplyr::summarise(
AGB_Mgha = sum(AGB_Mgha),
.groups = "drop") |>
dplyr::group_by(ano, Tratamento, plantio)|>
dplyr::summarise(meanAGB = mean(AGB_Mgha),
sd = sd(AGB_Mgha),
n = dplyr::n(),
se = sd/sqrt(n),
.groups = "drop") |>
dplyr::mutate(
Ano = as.factor(ifelse(ano == 2022, "5", "7")),
Tratamento = forcats::fct_recode(Tratamento,
"Nuclei\n25%" = "Nuclei 25%",
"Strips\n25%" = "Strips 25%",
"Nuclei\n50%" = "Nuclei 50%",
"Strips\n50%" = "Strips 50%",
"Planting\n100%" = "Planting 100%"
))
max_year <- agb_plot_data |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(meanAGB = max(meanAGB),.groups="drop")|>
dplyr::group_by(ano)|>
dplyr::summarise(max_year = max(meanAGB), .groups="drop")
agb_plot_data |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(meanAGB = max(meanAGB),.groups="drop")|>
dplyr::right_join(max_year, by="ano") |>
dplyr::mutate(porcentos = round( (meanAGB/max_year)*100, 0))# A tibble: 2 × 5
ano Tratamento meanAGB max_year porcentos
<dbl> <ord> <dbl> <dbl> <dbl>
1 2022 <NA> 48.9 48.9 100
2 2024 <NA> 83.2 83.2 100
Plotting
error_bars <- agb_plot_data |>
arrange(Ano, Tratamento, desc(plantio)) |>
group_by(Ano, Tratamento) |>
mutate(meanAGB_new = cumsum(meanAGB)) |>
ungroup()
porcentos <- data.frame(
Ano = rep(c("5", "7"), each = 5),
Tratamento = rep(unique(agb_plot_data$Tratamento), times = length(unique(agb_plot_data$Ano))),
meanAGB = c(8, 9, 12, 15, 27, 12, 17, 23, 25, 38),
label = c(c("17%","27%", "39%", "46%", "100%",
"21%", "34%", "52%", "53%", "100%"))
)
text_data <- data.frame(
Ano = rep(c("5", "7"), each = 5),
Tratamento = rep(unique(agb_plot_data$Tratamento), times = length(unique(agb_plot_data$Ano))),
meanAGB = c(10, 11, 14, 17, 29, 14, 19, 25, 27, 40),
label = c(c("a","a", "a", "a", "b",
"A", "AB", "B", "B", "C"))
)
fig_biomass<-agb_plot_data |>
ggplot(aes(x = Ano, y = meanAGB, fill = plantio)) +
geom_bar(stat = "identity", color = "black") +
geom_errorbar(data=error_bars,
aes(ymin = meanAGB_new - se, ymax = meanAGB_new + se),
width = .3) +
scale_fill_manual(
labels = c('Planted', 'Unplanted'),
values = c("#92D050", "white")
) +
labs(fill = "",
y = expression("Total aboveground biomass (Mg.ha"^{-1} ~")"),
x = "") +
facet_grid(.~Tratamento, switch = "both")+
geom_text(data = text_data, aes(x = Ano, y = meanAGB, label = label), size = 3, color = "black", inherit.aes = FALSE) +
geom_text(data = porcentos, aes(x = Ano, y = meanAGB, label = label), size = 3, color = "black", inherit.aes = FALSE) +
scale_y_continuous(expand = c(0, 0), limits = c(0,43))+
ggthemes::theme_clean()+
theme(
#plot.margin = margin(b = 0, l = 0),
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_blank(),
axis.text.x = element_text(color="black", hjust = 1, size=10),
axis.text.y = element_text(color="black", size=9),
axis.title.y = element_text(color="black", size=11),
strip.text = element_text(color="black", hjust = 0.5, size=10),
panel.spacing = unit(0.1, "mm"),
strip.placement = "outside",
legend.text = element_text(size = 11),
legend.title = element_blank(),
legend.key.size = unit(3, "mm"),
legend.box.background = element_blank(),
legend.box.margin = margin(t = 5, r = 0, b = 0, l = 0, "mm"),
legend.background = element_blank(),
legend.margin = margin(0,0,0,0),
legend.position = c(0.18, 0.90)
)
fig_biomassWarning: Removed 2 rows containing missing values or values outside the scale range
(`geom_bar()`).
# ggsave("fig_biomass.png", plot=fig_biomass, width=90, height = 90 ,unit="mm")Canopy Height Metric (CHM_mean, CHM_SD, and CHM_CV), Canopy cover (Canopy_cover), Rugosity (Rugosity_SD, and Rugosity_CV), Leaf Area Index (LAI_mean, LAI_SD, and LAI_CV), and LAIunder (LAIunder_mean, LAIunder_SD, LAIunder_CV).
2022
CHM2022 <- lidar_wide |>
dplyr::filter(ano==2022)|>
with(lme4::lmer(CHM_mean ~ Tratamento + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(CHM2022)
DHARMa::plotQQunif(CHM2022)2024
CHM2024 <- lidar_wide |>
dplyr::filter(ano==2024)|>
with(lme4::lmer(CHM_mean ~ Tratamento + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(CHM2024)
DHARMa::plotQQunif(CHM2024)Together
modeloCHM <- lidar_wide |>
with(lme4::lmer(CHM_mean ~ Tratamento*ano + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(modeloCHM)
DHARMa::plotQQunif(modeloCHM)Comparing estimated means
emmeans::emmeans(modeloCHM, pairwise ~ Tratamento | ano) |> multcomp::cld()ano = 2022:
Tratamento emmean SE df lower.CL upper.CL .group
Nuclei 25% 1.89 0.437 10 0.919 2.87 1
Strips 25% 2.12 0.437 10 1.142 3.09 1
Strips 50% 2.93 0.437 10 1.956 3.90 1
Nuclei 50% 3.00 0.437 10 2.027 3.97 1
Planting 100% 4.91 0.437 10 3.941 5.89 2
ano = 2024:
Tratamento emmean SE df lower.CL upper.CL .group
Nuclei 25% 2.76 0.437 10 1.788 3.73 1
Strips 25% 3.16 0.437 10 2.187 4.13 1
Strips 50% 4.43 0.437 10 3.454 5.40 2
Nuclei 50% 4.70 0.437 10 3.724 5.67 2
Planting 100% 6.99 0.437 10 6.018 7.96 3
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Comparing estimated slopes
emmeans::emtrends(modeloCHM, pairwise ~ Tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld() Tratamento ano.trend SE df lower.CL upper.CL .group
Nuclei 25% 0.434 0.201 36 0.026 0.843 1
Strips 25% 0.523 0.201 36 0.114 0.931 1
Strips 50% 0.749 0.201 36 0.341 1.157 1
Nuclei 50% 0.849 0.201 36 0.440 1.257 1
Planting 100% 1.039 0.201 36 0.630 1.447 1
Results are averaged over the levels of: ano
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Comparing fixed and random factors variance
MuMIn::r.squaredGLMM(modeloCHM) R2m R2c
[1,] 0.7002626 0.8728069
0.7002626 /0.8728069[1] 0.802311
Results:
The trend of variation around the mean CHM was the same over time among treatments (n.s.). The significant differences (\(p<0.001\)) observed between the 100% planting treatment and the other treatments in 2022 persisted in 2024.
sjPlot::tab_model(modeloCHM)| Dependent variable | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -876.29 | -1700.80 – -51.79 | 0.038 |
| Tratamento [Strips 25%] | -178.32 | -1344.35 – 987.70 | 0.759 |
| Tratamento [Nuclei 50%] | -836.96 | -2002.98 – 329.07 | 0.154 |
| Tratamento [Strips 50%] | -635.54 | -1801.57 – 530.48 | 0.277 |
| Tratamento [Planting 100%] |
-1218.82 | -2384.85 – -52.80 | 0.041 |
| ano | 0.43 | 0.03 – 0.84 | 0.037 |
| Tratamento [Strips 25%] × ano |
0.09 | -0.49 – 0.66 | 0.758 |
| Tratamento [Nuclei 50%] × ano |
0.41 | -0.16 – 0.99 | 0.154 |
| Tratamento [Strips 50%] × ano |
0.31 | -0.26 – 0.89 | 0.276 |
| Tratamento [Planting 100%] × ano |
0.60 | 0.03 – 1.18 | 0.040 |
| Random Effects | |||
| σ2 | 0.41 | ||
| τ00 Bloco | 0.55 | ||
| ICC | 0.58 | ||
| N Bloco | 5 | ||
| Observations | 50 | ||
| Marginal R2 / Conditional R2 | 0.700 / 0.873 | ||
broom.mixed::tidy(modeloCHM) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 12 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) -876. 407. -2.15
2 fixed <NA> TratamentoStrips 25% -178. 576. -0.31
3 fixed <NA> TratamentoNuclei 50% -837. 576. -1.45
4 fixed <NA> TratamentoStrips 50% -636. 576. -1.10
5 fixed <NA> TratamentoPlanting 100% -1219. 576. -2.12
6 fixed <NA> ano 0.434 0.201 2.16
7 fixed <NA> TratamentoStrips 25%:ano 0.088 0.285 0.31
8 fixed <NA> TratamentoNuclei 50%:ano 0.414 0.285 1.46
9 fixed <NA> TratamentoStrips 50%:ano 0.315 0.285 1.11
10 fixed <NA> TratamentoPlanting 100%:ano 0.604 0.285 2.12
11 ran_pars Bloco sd__(Intercept) 0.742 NA NA
12 ran_pars Residual sd__Observation 0.637 NA NA
Values per year
lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = round(mean(CHM_mean), 2),
SE = round(sd(CHM_mean)/sqrt(dplyr::n()), 2)) |>
dplyr::mutate(values = paste(Mean, "pm", SE)) |>
dplyr::select(-c(Mean, SE)) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)`summarise()` has grouped output by 'ano'. You can override using the `.groups`
argument.
# A tibble: 2 × 6
# Groups: ano [2]
ano `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<dbl> <chr> <chr> <chr> <chr> <chr>
1 2022 1.89 pm 0.3 2.12 pm 0.41 3 pm 0.26 2.93 pm 0.46 4.91 pm 0.33
2 2024 2.76 pm 0.48 3.16 pm 0.56 4.7 pm 0.28 4.43 pm 0.69 6.99 pm 0.4
“Slope” (increment/year)
lidar_wide |>
dplyr::group_by(Tratamento, ano) |>
dplyr::summarise(
media = mean(CHM_mean, na.rm = TRUE),
se = sd(CHM_mean, na.rm = TRUE) / sqrt(dplyr::n()),
.groups = "drop"
) |>
tidyr::pivot_wider(
names_from = ano,
values_from = c(media, se),
names_sep = "_"
) |>
dplyr::mutate(
variacao_media = round((media_2024 - media_2022) / 2, 2),
SE_variacao = round(sqrt(se_2022^2 + se_2024^2) / 2, 2)) |>
dplyr::mutate(values = paste(variacao_media, "pm", SE_variacao)) |>
dplyr::select(Tratamento, values) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)# A tibble: 1 × 5
`Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<chr> <chr> <chr> <chr> <chr>
1 0.43 pm 0.28 0.52 pm 0.35 0.85 pm 0.19 0.75 pm 0.42 1.04 pm 0.26
lidar_wide |>
dplyr::group_by(ano)|>
dplyr::summarise(mean(Canopy_cover),
min(Canopy_cover),
max(Canopy_cover),
sd(Canopy_cover))# A tibble: 2 × 5
ano `mean(Canopy_cover)` `min(Canopy_cover)` `max(Canopy_cover)`
<dbl> <dbl> <dbl> <dbl>
1 2022 38.4 5.29 99.8
2 2024 67.6 10.7 99.8
# ℹ 1 more variable: `sd(Canopy_cover)` <dbl>
2022
CC2022 <- lidar_wide |>
dplyr::filter(ano==2022)|>
with( glmmTMB::glmmTMB(Canopy_cover ~ Tratamento + (1|Bloco)) )Warning in glmmTMB::glmmTMB(Canopy_cover ~ Tratamento + (1 | Bloco)): use of
the 'data' argument is recommended
CC2022a <- lidar_wide |>
dplyr::filter(ano==2022)|>
with(glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento + (1 | Bloco),
family = glmmTMB::beta_family()))Warning in glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento + (1 | Bloco), : use
of the 'data' argument is recommended
AIC(CC2022, CC2022a) df AIC
CC2022 7 228.16951
CC2022a 7 -5.60719
par(mfrow=c(2,2))
DHARMa::plotResiduals(CC2022a)
DHARMa::plotQQunif(CC2022a)
DHARMa::plotResiduals(CC2022)
DHARMa::plotQQunif(CC2022)2024
CC2024 <- lidar_wide |>
dplyr::filter(ano==2024) |>
with(lme4::lmer(Canopy_cover ~ Tratamento + (1|Bloco) ) )
CC2024a <- lidar_wide |>
dplyr::filter(ano==2024) |>
with(glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento + (1 | Bloco),
family = glmmTMB::beta_family()))Warning in glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento + (1 | Bloco), : use
of the 'data' argument is recommended
AIC(CC2024, CC2024a) df AIC
CC2024 7 212.04388
CC2024a 7 -16.11127
par(mfrow=c(2,2))
DHARMa::plotResiduals(CC2024a)
DHARMa::plotQQunif(CC2024a)
DHARMa::plotResiduals(CC2024)
DHARMa::plotQQunif(CC2024)models
modeloCC <- lidar_wide |>
with(lme4::lmer(Canopy_cover ~ Tratamento*ano + (1|Bloco) ) )
modeloCCa <- lidar_wide |>
with(glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento*ano + (1 | Bloco),
family = glmmTMB::beta_family()))Warning in glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento * ano + (1 | Bloco),
: use of the 'data' argument is recommended
Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
NA/NaN function evaluation
Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
NA/NaN function evaluation
par(mfrow=c(2,2))
DHARMa::plotResiduals(modeloCC)
DHARMa::plotQQunif(modeloCC)
DHARMa::plotResiduals(modeloCCa)
DHARMa::plotQQunif(modeloCCa)1 outlier
Estimated means
emmeans::emmeans(modeloCC, pairwise ~ Tratamento | ano) |> multcomp::cld()ano = 2022:
Tratamento emmean SE df lower.CL upper.CL .group
Strips 25% 23.7 11.5 20.8 -0.25 47.6 1
Nuclei 50% 26.5 11.5 20.8 2.56 50.4 1
Strips 50% 31.3 11.5 20.8 7.42 55.3 1
Nuclei 25% 31.6 11.5 20.8 7.70 55.5 1
Planting 100% 79.1 11.5 20.8 55.14 103.0 2
ano = 2024:
Tratamento emmean SE df lower.CL upper.CL .group
Strips 25% 44.7 11.5 20.8 20.78 68.6 1
Nuclei 25% 56.5 11.5 20.8 32.62 80.5 1
Strips 50% 69.8 11.5 20.8 45.88 93.7 12
Nuclei 50% 70.9 11.5 20.8 47.02 94.9 12
Planting 100% 96.2 11.5 20.8 72.29 120.1 2
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Estimated slopes
emmeans::emtrends(modeloCC, pairwise ~ Tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld() Tratamento ano.trend SE df lower.CL upper.CL .group
Planting 100% 8.57 6.71 36 -5.02 22.2 1
Strips 25% 10.52 6.71 36 -3.08 24.1 1
Nuclei 25% 12.46 6.71 36 -1.14 26.1 1
Strips 50% 19.23 6.71 36 5.63 32.8 1
Nuclei 50% 22.23 6.71 36 8.63 35.8 1
Results are averaged over the levels of: ano
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Comparing random and fixed effects
MuMIn::r.squaredGLMM(modeloCC) R2m R2c
[1,] 0.469023 0.6387404
0.469023 / 0.6387404[1] 0.7342936
sjPlot::tab_model(modeloCC)| Dependent variable | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -25156.96 | -52617.32 – 2303.39 | 0.071 |
| Tratamento [Strips 25%] | 3917.03 | -34917.77 – 42751.83 | 0.839 |
| Tratamento [Nuclei 50%] | -19762.19 | -58596.99 – 19072.61 | 0.309 |
| Tratamento [Strips 50%] | -13690.08 | -52524.88 – 25144.72 | 0.480 |
| Tratamento [Planting 100%] |
7898.75 | -30936.05 – 46733.55 | 0.683 |
| ano | 12.46 | -1.12 – 26.03 | 0.071 |
| Tratamento [Strips 25%] × ano |
-1.94 | -21.14 – 17.26 | 0.839 |
| Tratamento [Nuclei 50%] × ano |
9.77 | -9.43 – 28.97 | 0.309 |
| Tratamento [Strips 50%] × ano |
6.77 | -12.43 – 25.97 | 0.480 |
| Tratamento [Planting 100%] × ano |
-3.88 | -23.08 – 15.31 | 0.684 |
| Random Effects | |||
| σ2 | 449.60 | ||
| τ00 Bloco | 211.22 | ||
| ICC | 0.32 | ||
| N Bloco | 5 | ||
| Observations | 50 | ||
| Marginal R2 / Conditional R2 | 0.469 / 0.639 | ||
broom.mixed::tidy(modeloCC) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 12 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) -25157. 13565. -1.86
2 fixed <NA> TratamentoStrips 25% 3917. 19183. 0.204
3 fixed <NA> TratamentoNuclei 50% -19762. 19183. -1.03
4 fixed <NA> TratamentoStrips 50% -13690. 19183. -0.714
5 fixed <NA> TratamentoPlanting 100% 7899. 19183. 0.412
6 fixed <NA> ano 12.5 6.70 1.86
7 fixed <NA> TratamentoStrips 25%:ano -1.94 9.48 -0.205
8 fixed <NA> TratamentoNuclei 50%:ano 9.77 9.48 1.03
9 fixed <NA> TratamentoStrips 50%:ano 6.77 9.48 0.714
10 fixed <NA> TratamentoPlanting 100%:ano -3.88 9.48 -0.409
11 ran_pars Bloco sd__(Intercept) 14.5 NA NA
12 ran_pars Residual sd__Observation 21.2 NA NA
per year
lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = round(mean(Canopy_cover), 2),
SE = round((sd(Canopy_cover)/sqrt(dplyr::n())), 2),
.groups = "drop") |>
dplyr::mutate(values = paste(Mean, "pm", SE)) |>
dplyr::select(-c(Mean, SE)) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)# A tibble: 2 × 6
ano `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<dbl> <chr> <chr> <chr> <chr> <chr>
1 2022 31.62 pm 17.23 23.67 pm 4.99 26.48 pm 1.81 31.34 pm 7.74 79.06 pm 10.54
2 2024 56.54 pm 18.23 44.7 pm 13.14 70.94 pm 11.3 69.8 pm 13.76 96.2 pm 1.89
“Slope” (increment/year)
lidar_wide |>
dplyr::group_by(Tratamento, ano) |>
dplyr::summarise(
media = mean(Canopy_cover, na.rm = TRUE),
se = sd(Canopy_cover, na.rm = TRUE) / sqrt(dplyr::n()),
.groups = "drop"
) |>
tidyr::pivot_wider(
names_from = ano,
values_from = c(media, se),
names_sep = "_"
) |>
dplyr::mutate(
variacao_media = round((media_2024 - media_2022) / 2, 2),
SE_variacao = round(sqrt(se_2022^2 + se_2024^2) / 2, 2)) |>
dplyr::mutate(values = paste(variacao_media, "pm", SE_variacao)) |>
dplyr::select(Tratamento, values) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)# A tibble: 1 × 5
`Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<chr> <chr> <chr> <chr> <chr>
1 12.46 pm 12.54 10.52 pm 7.03 22.23 pm 5.72 19.23 pm 7.89 8.57 pm 5.35
Performing a pairwise comparison of variances is not as straightforward as comparing means, but you can use specific statistical techniques to compare variances between treatment groups. Here’s how you can do this in R:
pairwise_levene <- function(data, response, group) {
treatment_pairs <- combn(unique(data[[group]]), 2, simplify = FALSE)
p_values <- vector("list", length(treatment_pairs))
for (i in seq_along(treatment_pairs)) { #i =1
subset_data <- data[data[[group]] %in% treatment_pairs[[i]], ]
test_result <- leveneTest(subset_data[[response]] ~ subset_data[[group]])
p_values[[i]] <- test_result[1, "Pr(>F)"]
}
comparison <- sapply(treatment_pairs, function(x) paste(x, collapse = " vs "))
results <- data.frame(Comparison = comparison, P_Value = unlist(p_values))
results$Adjusted_P_Value <- p.adjust(results$P_Value, method = "bonferroni")
return(results)
}
pairwise_levene(lidar2022_wide , "Canopy_cover", "Tratamento")Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Comparison P_Value Adjusted_P_Value
1 25_faixa vs 25_tabuleiro 0.41911844 1.0000000
2 25_faixa vs 50_faixa 0.60465350 1.0000000
3 25_faixa vs 50_tabuleiro 0.52704716 1.0000000
4 25_faixa vs 100_controle 0.21179066 1.0000000
5 25_tabuleiro vs 50_faixa 0.56742080 1.0000000
6 25_tabuleiro vs 50_tabuleiro 0.31075886 1.0000000
7 25_tabuleiro vs 100_controle 0.89320378 1.0000000
8 50_faixa vs 50_tabuleiro 0.26004211 1.0000000
9 50_faixa vs 100_controle 0.43307531 1.0000000
10 50_tabuleiro vs 100_controle 0.07724271 0.7724271
pairwise_levene(lidar2024_wide , "Canopy_cover", "Tratamento")Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Comparison P_Value Adjusted_P_Value
1 25_faixa vs 25_tabuleiro 0.47523130 1.0000000
2 25_faixa vs 50_faixa 0.73908327 1.0000000
3 25_faixa vs 50_tabuleiro 0.87724043 1.0000000
4 25_faixa vs 100_controle 0.28751774 1.0000000
5 25_tabuleiro vs 50_faixa 0.65306786 1.0000000
6 25_tabuleiro vs 50_tabuleiro 0.49847386 1.0000000
7 25_tabuleiro vs 100_controle 0.08104640 0.8104640
8 50_faixa vs 50_tabuleiro 0.81722803 1.0000000
9 50_faixa vs 100_controle 0.11052807 1.0000000
10 50_tabuleiro vs 100_controle 0.09110862 0.9110862
lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = round(mean(Rugosity_SD),2),
SE = round((sd(Rugosity_SD)/sqrt(dplyr::n())), 2),
.groups="drop") |>
dplyr::mutate(values = paste(Mean, "pm", SE)) |>
dplyr::select(-c(Mean, SE)) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)# A tibble: 2 × 6
ano `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<dbl> <chr> <chr> <chr> <chr> <chr>
1 2022 2.08 pm 0.22 2.38 pm 0.17 2.73 pm 0.31 2.64 pm 0.19 2.57 pm 0.06
2 2024 2.77 pm 0.32 3.15 pm 0.18 3.41 pm 0.36 3.3 pm 0.17 2.79 pm 0.09
Increment/year
lidar_wide |>
dplyr::group_by(Tratamento, ano) |>
dplyr::summarise(
media = mean(Rugosity_SD, na.rm = TRUE),
se = sd(Rugosity_SD, na.rm = TRUE) / sqrt(dplyr::n()),
.groups = "drop"
) |>
tidyr::pivot_wider(
names_from = ano,
values_from = c(media, se),
names_sep = "_"
) |>
dplyr::mutate(
variacao_media = round((media_2024 - media_2022) / 2, 2),
SE_variacao = round(sqrt(se_2022^2 + se_2024^2) / 2, 2)) |>
dplyr::mutate(values = paste(variacao_media, "pm", SE_variacao)) |>
dplyr::select(Tratamento, values) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)# A tibble: 1 × 5
`Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<chr> <chr> <chr> <chr> <chr>
1 0.34 pm 0.19 0.39 pm 0.12 0.34 pm 0.24 0.33 pm 0.13 0.11 pm 0.05
2022
LAI2022 <- lidar_wide |>
dplyr::filter(ano==2022) |>
with(lme4::lmer(LAI_mean ~ Tratamento + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(LAI2022)
DHARMa::plotQQunif(LAI2022)2024
LAI2024 <- lidar_wide |>
dplyr::filter(ano==2024) |>
with(lme4::lmer(LAI_mean ~ Tratamento + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(LAI2024)
DHARMa::plotQQunif(LAI2024)Full model
modeloLAI <- lidar_wide |>
with(lme4::lmer(LAI_mean ~ Tratamento * ano + (1 | Bloco)))
par(mfrow=c(1,2))
DHARMa::plotResiduals(modeloLAI)
DHARMa::plotQQunif(modeloLAI)Pairwise comparison
emmeans::emmeans(modeloLAI, pairwise ~ Tratamento | ano) |> multcomp::cld()ano = 2022:
Tratamento emmean SE df lower.CL upper.CL .group
Nuclei 25% 0.447 0.159 13.2 0.105 0.789 1
Strips 25% 0.591 0.159 13.2 0.249 0.933 1
Nuclei 50% 0.781 0.159 13.2 0.439 1.123 1
Strips 50% 0.828 0.159 13.2 0.485 1.170 1
Planting 100% 1.565 0.159 13.2 1.223 1.907 2
ano = 2024:
Tratamento emmean SE df lower.CL upper.CL .group
Nuclei 25% 0.900 0.159 13.2 0.558 1.242 1
Strips 25% 1.047 0.159 13.2 0.705 1.389 12
Strips 50% 1.353 0.159 13.2 1.011 1.695 12
Nuclei 50% 1.383 0.159 13.2 1.041 1.725 2
Planting 100% 2.408 0.159 13.2 2.066 2.750 3
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Trend:
emmeans::emtrends(modeloLAI, pairwise ~ Tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld() Tratamento ano.trend SE df lower.CL upper.CL .group
Nuclei 25% 0.226 0.0812 36 0.0617 0.391 1
Strips 25% 0.228 0.0812 36 0.0633 0.393 1
Strips 50% 0.263 0.0812 36 0.0980 0.428 1
Nuclei 50% 0.301 0.0812 36 0.1360 0.466 1
Planting 100% 0.421 0.0812 36 0.2565 0.586 1
Results are averaged over the levels of: ano
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
Fixed vs random
MuMIn::r.squaredGLMM(modeloLAI) R2m R2c
[1,] 0.7064328 0.8459851
0.7064328 /0.8459851[1] 0.8350417
sjPlot::tab_model(modeloLAI)| Dependent variable | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -457.36 | -790.02 – -124.71 | 0.008 |
| Tratamento [Strips 25%] | -3.08 | -473.52 – 467.36 | 0.989 |
| Tratamento [Nuclei 50%] | -150.02 | -620.46 – 320.42 | 0.522 |
| Tratamento [Strips 50%] | -73.12 | -543.56 – 397.32 | 0.755 |
| Tratamento [Planting 100%] |
-392.88 | -863.32 – 77.56 | 0.099 |
| ano | 0.23 | 0.06 – 0.39 | 0.008 |
| Tratamento [Strips 25%] × ano |
0.00 | -0.23 – 0.23 | 0.989 |
| Tratamento [Nuclei 50%] × ano |
0.07 | -0.16 – 0.31 | 0.521 |
| Tratamento [Strips 50%] × ano |
0.04 | -0.20 – 0.27 | 0.753 |
| Tratamento [Planting 100%] × ano |
0.19 | -0.04 – 0.43 | 0.098 |
| Random Effects | |||
| σ2 | 0.07 | ||
| τ00 Bloco | 0.06 | ||
| ICC | 0.48 | ||
| N Bloco | 5 | ||
| Observations | 50 | ||
| Marginal R2 / Conditional R2 | 0.706 / 0.846 | ||
broom.mixed::tidy(modeloLAI) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 12 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) -457. 164. -2.78
2 fixed <NA> TratamentoStrips 25% -3.08 232. -0.013
3 fixed <NA> TratamentoNuclei 50% -150. 232. -0.646
4 fixed <NA> TratamentoStrips 50% -73.1 232. -0.315
5 fixed <NA> TratamentoPlanting 100% -393. 232. -1.69
6 fixed <NA> ano 0.226 0.081 2.79
7 fixed <NA> TratamentoStrips 25%:ano 0.002 0.115 0.014
8 fixed <NA> TratamentoNuclei 50%:ano 0.074 0.115 0.647
9 fixed <NA> TratamentoStrips 50%:ano 0.036 0.115 0.316
10 fixed <NA> TratamentoPlanting 100%:ano 0.195 0.115 1.70
11 ran_pars Bloco sd__(Intercept) 0.245 NA NA
12 ran_pars Residual sd__Observation 0.257 NA NA
lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = round(mean(LAI_mean),2),
SE = round((sd(LAI_mean)/sqrt(dplyr::n())), 2),
.groups = "drop") |>
dplyr::mutate(values = paste(Mean, "pm", SE)) |>
dplyr::select(-c(Mean, SE)) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)# A tibble: 2 × 6
ano `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<dbl> <chr> <chr> <chr> <chr> <chr>
1 2022 0.45 pm 0.09 0.59 pm 0.12 0.78 pm 0.08 0.83 pm 0.16 1.57 pm 0.17
2 2024 0.9 pm 0.15 1.05 pm 0.2 1.38 pm 0.16 1.35 pm 0.19 2.41 pm 0.21
Increment/year (slope)
lidar_wide |>
dplyr::group_by(Tratamento, ano) |>
dplyr::summarise(
media = mean(LAI_mean, na.rm = TRUE),
se = sd(LAI_mean, na.rm = TRUE) / sqrt(dplyr::n()),
.groups = "drop"
) |>
tidyr::pivot_wider(
names_from = ano,
values_from = c(media, se),
names_sep = "_"
) |>
dplyr::mutate(
variacao_media = round((media_2024 - media_2022) / 2, 2),
SE_variacao = round(sqrt(se_2022^2 + se_2024^2) / 2, 2)) |>
dplyr::mutate(values = paste(variacao_media, "pm", SE_variacao)) |>
dplyr::select(Tratamento, values) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)# A tibble: 1 × 5
`Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<chr> <chr> <chr> <chr> <chr>
1 0.23 pm 0.09 0.23 pm 0.12 0.3 pm 0.09 0.26 pm 0.12 0.42 pm 0.14
max_year <- lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(CHM_mean),.groups="drop") |>
dplyr::group_by(ano)|>
dplyr::summarise(max_year = max(Mean),.groups="drop")
lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(CHM_mean),.groups="drop") |>
dplyr::right_join(max_year, by="ano") |>
dplyr::mutate(porcentos = round( (Mean/max_year)*100, 0))# A tibble: 10 × 5
ano Tratamento Mean max_year porcentos
<dbl> <fct> <dbl> <dbl> <dbl>
1 2022 Nuclei 25% 1.89 4.91 39
2 2022 Strips 25% 2.12 4.91 43
3 2022 Nuclei 50% 3.00 4.91 61
4 2022 Strips 50% 2.93 4.91 60
5 2022 Planting 100% 4.91 4.91 100
6 2024 Nuclei 25% 2.76 6.99 39
7 2024 Strips 25% 3.16 6.99 45
8 2024 Nuclei 50% 4.70 6.99 67
9 2024 Strips 50% 4.43 6.99 63
10 2024 Planting 100% 6.99 6.99 100
Plot
letras <- data.frame(
Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
ano = c(2022, 2024),
label = c("A "," a",
"A "," a",
"A "," b",
"A "," b",
"B "," c"),
y = c(2.8, 3.9,
3.1, 4.5,
3.9, 5.7,
3.9, 5.7,
5.9, 7.8) )
porcentos <- data.frame(
Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
ano = c(2022, 2024),
label = c("39% "," 39%",
"43% "," 45%",
"61% "," 67%",
"60% "," 63%",
"100% "," 100%"),
y = c(2.2, 3.1,
2.5, 3.7,
3.3, 5,
3.3, 5,
5.3, 7.1)
)
p1 <- lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(CHM_mean),
SD = sd(CHM_mean),
SE = SD/sqrt(dplyr::n()),
.groups="drop") |>
ggplot(aes(x = Tratamento, y = Mean, fill = as.factor(ano))) +
geom_bar(position = position_dodge2(preserve = "single"),
stat = "identity", color = "black") +
geom_errorbar(aes(ymin = Mean - SE, ymax = Mean + SE), width = .2, position = position_dodge(0.9)) +
scale_fill_manual(
labels = c('5 years', '7 years'),
values = c("#CAF0F8", "#00B4D8")) +
geom_text(data = letras, aes(x = Tratamento, y = y, label = label),size = 3.5, color="black", vjust = -0.5) +
geom_text(data = porcentos, aes(x = Tratamento, y = y, label = label),size = 3, color="black", vjust = -0.5) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 9), breaks = seq(0, 8, by = 1))+
labs(fill = "",
y = "Canopy height model (m)",
x = "") +
theme_classic()+
theme(
axis.text.x = element_blank(),
axis.text.y = element_text(color="black", size=9),
axis.title.y = element_text(color="black", size=10),
legend.position = c(0.15, 0.9),
legend.text = element_text(size = 10),
legend.key.size = unit(3, "mm"),
plot.margin = margin(t=5, 0, 0, 0, "mm"))
p1max_year <- lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(Canopy_cover),.groups="drop") |>
dplyr::group_by(ano)|>
dplyr::summarise(max_year = max(Mean),.groups="drop")
lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(Canopy_cover),.groups="drop") |>
dplyr::right_join(max_year, by="ano") |>
dplyr::mutate(porcentos = round( (Mean/max_year)*100, 0))# A tibble: 10 × 5
ano Tratamento Mean max_year porcentos
<dbl> <fct> <dbl> <dbl> <dbl>
1 2022 Nuclei 25% 31.6 79.1 40
2 2022 Strips 25% 23.7 79.1 30
3 2022 Nuclei 50% 26.5 79.1 33
4 2022 Strips 50% 31.3 79.1 40
5 2022 Planting 100% 79.1 79.1 100
6 2024 Nuclei 25% 56.5 96.2 59
7 2024 Strips 25% 44.7 96.2 46
8 2024 Nuclei 50% 70.9 96.2 74
9 2024 Strips 50% 69.8 96.2 73
10 2024 Planting 100% 96.2 96.2 100
plot
letras <- data.frame(
Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
ano = c(2022, 2024),
label = c("A "," a",
"A "," a",
"A "," ab",
"A "," ab",
"B "," b"),
y = c(52, 77,
40, 65,
40, 85,
45, 90,
97, 109)
)
porcentos <- data.frame(
Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
ano = c(2022, 2024),
label = c("40% "," 59%",
"30% "," 46%",
"33% ", " 74%",
"40% "," 73%",
"100% "," 100%"),
y = c(42, 67,
30, 55,
30, 75,
35, 80,
87, 100)
)
p2 <- lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(Canopy_cover),
SD = sd(Canopy_cover),
SE = SD/sqrt(dplyr::n()),
.groups="drop") |>
ggplot(aes(x = Tratamento, y = Mean, fill = as.factor(ano))) +
geom_bar(position = position_dodge2(preserve = "single"),
stat = "identity", color = "black") +
geom_errorbar(aes(ymin = Mean - SE, ymax = Mean + SE), width = .2, position = position_dodge(0.9)) +
scale_fill_manual(
labels = c('5 years', '7 years'),
values = c("#CAF0F8", "#00B4D8")) +
geom_text(data = letras, aes(x = Tratamento, y = y, label = label),size = 3.5, color="black", vjust = -0.5) +
geom_text(data = porcentos, aes(x = Tratamento, y = y, label = label),size = 3, color="black", vjust = -0.5) +
scale_y_continuous(expand = c(0, 0), limits=c(0, 120), breaks = seq(0, 100, by = 20))+
labs(fill = "",
y = "Canopy cover (%)",
x = "") +
theme_classic()+
theme(
axis.text.x = element_blank(),
axis.text.y = element_text(color="black", size=10),
axis.title.y = element_text(color="black", size=10),
plot.margin = margin(t=5, 0, 0, 0, "mm"),
legend.position = 'none')
p2max_year <- lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(LAI_mean),.groups="drop") |>
dplyr::group_by(ano)|>
dplyr::summarise(max_year = max(Mean),.groups="drop")
lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(LAI_mean),.groups="drop") |>
dplyr::right_join(max_year, by="ano") |>
dplyr::mutate(porcentos = round( (Mean/max_year)*100, 0))# A tibble: 10 × 5
ano Tratamento Mean max_year porcentos
<dbl> <fct> <dbl> <dbl> <dbl>
1 2022 Nuclei 25% 0.447 1.57 29
2 2022 Strips 25% 0.591 1.57 38
3 2022 Nuclei 50% 0.781 1.57 50
4 2022 Strips 50% 0.828 1.57 53
5 2022 Planting 100% 1.57 1.57 100
6 2024 Nuclei 25% 0.900 2.41 37
7 2024 Strips 25% 1.05 2.41 43
8 2024 Nuclei 50% 1.38 2.41 57
9 2024 Strips 50% 1.35 2.41 56
10 2024 Planting 100% 2.41 2.41 100
Plot
letras <- data.frame(
Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
ano = c(2022, 2024),
label = c("A "," a",
"A "," ab",
"A ", " b",
"A "," ab",
"B "," c"),
y = c( 0.85, 1.4,
1.05, 1.55,
1.15, 1.82,
1.25, 1.82,
2.05, 2.9)
) |>
dplyr::mutate(
Tratamento = forcats::fct_recode(Tratamento,
"Nuclei\n25%" = "Nuclei 25%",
"Strips\n25%" = "Strips 25%",
"Nuclei\n50%" = "Nuclei 50%",
"Strips\n50%" = "Strips 50%",
"Planting\n100%" = "Planting 100%"
)
)
porcentos <- data.frame(
Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
ano = c(2022, 2024),
label = c("29% "," 37%",
"38% "," 43%",
"50% "," 57%",
"53% "," 56%",
"100% "," 100%"),
y = c( 0.55, 1.1,
0.75, 1.25,
0.85, 1.52,
0.95, 1.52,
1.75, 2.65)
)|>
dplyr::mutate(
Tratamento = forcats::fct_recode(Tratamento,
"Nuclei\n25%" = "Nuclei 25%",
"Strips\n25%" = "Strips 25%",
"Nuclei\n50%" = "Nuclei 50%",
"Strips\n50%" = "Strips 50%",
"Planting\n100%" = "Planting 100%"
)
)
p3 <- lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(LAI_mean),
SD = sd(LAI_mean),
SE = SD/sqrt(dplyr::n()),
.groups="drop") |>
dplyr::mutate(
Tratamento = forcats::fct_recode(Tratamento,
"Nuclei\n25%" = "Nuclei 25%",
"Strips\n25%" = "Strips 25%",
"Nuclei\n50%" = "Nuclei 50%",
"Strips\n50%" = "Strips 50%",
"Planting\n100%" = "Planting 100%"))|>
ggplot(aes(x = Tratamento, y = Mean, fill = as.factor(ano))) +
geom_bar(position = position_dodge2(preserve = "single"),
stat = "identity", color = "black") +
geom_errorbar(aes(ymin = Mean - SE, ymax = Mean + SE), width = .2, position = position_dodge(0.9)) +
scale_fill_manual(
labels = c('5 years', '7 years'),
values = c("#CAF0F8", "#00B4D8")) +
geom_text(data = letras, aes(x = Tratamento, y = y, label = label),size = 3.5, color="black", vjust = -0.5) +
geom_text(data = porcentos, aes(x = Tratamento, y = y, label = label),size = 3, color="black", vjust = -0.5) +
scale_y_continuous(expand = c(0, 0), limits = c(0,3.2), breaks = seq(0, 3, by=1))+
labs(fill = "",
y = expression("LAI (m"^{2} ~ "/ m"^{2}~ ")"),
x = "") +
theme_classic()+
theme(
axis.text.x = element_text(color="black", hjust = 0.5, size=10),
axis.text.y = element_text(color="black", size=10),
axis.title.y = element_text(color="black", size=10),
plot.margin = margin(t=5, 0, 0, 0, "mm"),
legend.position = 'none')
p3# 90 mm x 180 mm
fig_lidar <- cowplot::plot_grid(p1,p2,p3, ncol = 1,align = "v", axis = "l",
labels = c("(a)", "(b)", "(c)"), label_size = 10, label_x = 0.01, label_y = 1)
fig_lidar# ggsave("fig_lidar.png", plot=fig_lidar, width=90, height = 210 ,unit="mm")capital letters: 5 years
small letters: 7 years
Table of Costs (in R$)
costs <- data.frame(
tratamento = c("Planting 100%", "Strips 50%", "Nuclei 50%", "Strips 25%", "Nuclei 25%"),
labor_machinery = c(7860.00, 6407.90, 7563.90, 5131.05, 6395.23),
inputs = c(7032.50, 4096.25, 4096.25, 2628.13, 2628.13)) |>
dplyr::mutate(total = labor_machinery + inputs)cover_raw <- lidar_wide |>
dplyr::select(ano, Bloco, Tratamento, Canopy_cover) |>
dplyr::rename(bloco = Bloco, tratamento = Tratamento) |>
dplyr::mutate(tratamento = ordered(tratamento,
levels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%")),
bloco = ordered(bloco, levels = c("BI","BII", "BIII", "BIV", "BV")),
bloco = as.character(as.numeric(bloco)))|>
dplyr::left_join(costs) |>
dplyr::mutate(totalUS = (total/5)/5.8)Joining with `by = join_by(tratamento)`
cost_raw <- AGB_data |>
dplyr::select(ano, bloco, tratamento, AGB_Mgha) |>
dplyr::mutate(
tratamento = dplyr::case_when(
tratamento == "100" ~ "Planting 100%",
tratamento == "F25" ~ "Strips 25%",
tratamento == "F50" ~ "Strips 50%",
tratamento == "N25" ~ "Nuclei 25%",
tratamento == "N50" ~ "Nuclei 50%"
)
) |>
dplyr::right_join(cover_raw) |>
dplyr::mutate(agb_cost = totalUS/AGB_Mgha,
cc_cost = totalUS/Canopy_cover,
ano = ano-2017) |>
dplyr::select(ano:tratamento, agb_cost, cc_cost)Joining with `by = join_by(ano, bloco, tratamento)`
modeloCC_cost <- lme4::lmer(log(cc_cost) ~ tratamento*ano + (1|bloco), data=cost_raw)
par(mfrow=c(1,2))
DHARMa::plotResiduals(modeloCC_cost)
DHARMa::plotQQunif(modeloCC_cost)emmeans::emmeans(modeloCC_cost, pairwise ~ tratamento | ano) |> multcomp::cld()ano = 5:
tratamento emmean SE df lower.CL upper.CL .group
Planting 100% 1.91 0.256 18 1.37 2.45 1
Strips 25% 2.50 0.256 18 1.96 3.03 12
Strips 50% 2.56 0.256 18 2.02 3.10 12
Nuclei 50% 2.73 0.256 18 2.19 3.27 12
Nuclei 25% 2.77 0.256 18 2.24 3.31 2
ano = 7:
tratamento emmean SE df lower.CL upper.CL .group
Planting 100% 1.68 0.256 18 1.14 2.21 1
Strips 50% 1.75 0.256 18 1.21 2.29 1
Nuclei 50% 1.79 0.256 18 1.25 2.32 1
Strips 25% 1.92 0.256 18 1.38 2.46 1
Nuclei 25% 1.99 0.256 18 1.45 2.53 1
Degrees-of-freedom method: kenward-roger
Results are given on the log (not the response) scale.
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
sjPlot::tab_model(modeloCC_cost)| log(cc cost) | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 4.73 | 2.93 – 6.54 | <0.001 |
| tratamento [Nuclei 50%] | 0.36 | -2.15 – 2.86 | 0.775 |
| tratamento [Planting 100%] |
-2.23 | -4.74 – 0.27 | 0.079 |
| tratamento [Strips 25%] | -0.81 | -3.31 – 1.70 | 0.519 |
| tratamento [Strips 50%] | -0.15 | -2.66 – 2.36 | 0.905 |
| ano | -0.39 | -0.68 – -0.10 | 0.010 |
| tratamento [Nuclei 50%] × ano |
-0.08 | -0.49 – 0.33 | 0.697 |
| tratamento [Planting 100%] × ano |
0.27 | -0.14 – 0.69 | 0.186 |
| tratamento [Strips 25%] × ano |
0.11 | -0.31 – 0.52 | 0.607 |
| tratamento [Strips 50%] × ano |
-0.01 | -0.42 – 0.40 | 0.952 |
| Random Effects | |||
| σ2 | 0.21 | ||
| τ00 bloco | 0.12 | ||
| ICC | 0.37 | ||
| N bloco | 5 | ||
| Observations | 50 | ||
| Marginal R2 / Conditional R2 | 0.341 / 0.583 | ||
broom.mixed::tidy(modeloCC_cost) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), round, 3)) # A tibble: 12 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 4.73 0.89 5.32
2 fixed <NA> tratamentoNuclei 50% 0.356 1.24 0.287
3 fixed <NA> tratamentoPlanting 100% -2.23 1.24 -1.80
4 fixed <NA> tratamentoStrips 25% -0.806 1.24 -0.651
5 fixed <NA> tratamentoStrips 50% -0.149 1.24 -0.12
6 fixed <NA> ano -0.392 0.144 -2.72
7 fixed <NA> tratamentoNuclei 50%:ano -0.08 0.204 -0.392
8 fixed <NA> tratamentoPlanting 100%:ano 0.274 0.204 1.35
9 fixed <NA> tratamentoStrips 25%:ano 0.106 0.204 0.518
10 fixed <NA> tratamentoStrips 50%:ano -0.012 0.204 -0.061
11 ran_pars bloco sd__(Intercept) 0.347 NA NA
12 ran_pars Residual sd__Observation 0.455 NA NA
modeloAGB_cost <- lme4::lmer(log(agb_cost) ~ tratamento*ano + (1|bloco), data=cost_raw )
par(mfrow=c(1,2))
DHARMa::plotResiduals(modeloAGB_cost)
DHARMa::plotQQunif(modeloAGB_cost)emmeans::emmeans(modeloAGB_cost, pairwise ~ tratamento | ano) |> multcomp::cld()ano = 5:
tratamento emmean SE df lower.CL upper.CL .group
Planting 100% 3.21 0.217 17.5 2.75 3.67 1
Strips 50% 3.68 0.217 17.5 3.22 4.14 1
Strips 25% 3.76 0.217 17.5 3.30 4.22 1
Nuclei 50% 3.81 0.217 17.5 3.36 4.27 12
Nuclei 25% 4.50 0.217 17.5 4.04 4.95 2
ano = 7:
tratamento emmean SE df lower.CL upper.CL .group
Planting 100% 2.79 0.217 17.5 2.33 3.25 1
Strips 25% 3.08 0.217 17.5 2.63 3.54 1
Strips 50% 3.12 0.217 17.5 2.67 3.58 12
Nuclei 50% 3.17 0.217 17.5 2.71 3.62 12
Nuclei 25% 3.80 0.217 17.5 3.35 4.26 2
Degrees-of-freedom method: kenward-roger
Results are given on the log (not the response) scale.
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
And the slopes?
emmeans::emtrends(modeloCC_cost, pairwise ~ tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld() tratamento ano.trend SE df lower.CL upper.CL .group
Nuclei 50% -0.472 0.144 36 -0.764 -0.18000 1
Strips 50% -0.405 0.144 36 -0.697 -0.11250 1
Nuclei 25% -0.392 0.144 36 -0.684 -0.10005 1
Strips 25% -0.287 0.144 36 -0.579 0.00551 1
Planting 100% -0.118 0.144 36 -0.410 0.17419 1
Results are averaged over the levels of: ano
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
emmeans::emtrends(modeloAGB_cost, pairwise ~ tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld() tratamento ano.trend SE df lower.CL upper.CL .group
Nuclei 25% -0.346 0.121 36 -0.592 -0.1006 1
Strips 25% -0.338 0.121 36 -0.584 -0.0928 1
Nuclei 50% -0.324 0.121 36 -0.570 -0.0782 1
Strips 50% -0.280 0.121 36 -0.525 -0.0340 1
Planting 100% -0.210 0.121 36 -0.456 0.0355 1
Results are averaged over the levels of: ano
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
sjPlot::tab_model(modeloAGB_cost)| log(agb cost) | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 6.23 | 4.71 – 7.74 | <0.001 |
| tratamento [Nuclei 50%] | -0.79 | -2.90 – 1.32 | 0.452 |
| tratamento [Planting 100%] |
-1.96 | -4.07 – 0.15 | 0.067 |
| tratamento [Strips 25%] | -0.77 | -2.88 – 1.34 | 0.463 |
| tratamento [Strips 50%] | -1.15 | -3.26 – 0.96 | 0.278 |
| ano | -0.35 | -0.59 – -0.10 | 0.007 |
| tratamento [Nuclei 50%] × ano |
0.02 | -0.32 – 0.37 | 0.897 |
| tratamento [Planting 100%] × ano |
0.14 | -0.21 – 0.48 | 0.432 |
| tratamento [Strips 25%] × ano |
0.01 | -0.34 – 0.35 | 0.964 |
| tratamento [Strips 50%] × ano |
0.07 | -0.28 – 0.41 | 0.700 |
| Random Effects | |||
| σ2 | 0.15 | ||
| τ00 bloco | 0.09 | ||
| ICC | 0.38 | ||
| N bloco | 5 | ||
| Observations | 50 | ||
| Marginal R2 / Conditional R2 | 0.498 / 0.688 | ||
broom.mixed::tidy(modeloAGB_cost) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), round, 3)) # A tibble: 12 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 6.23 0.749 8.31
2 fixed <NA> tratamentoNuclei 50% -0.792 1.04 -0.76
3 fixed <NA> tratamentoPlanting 100% -1.96 1.04 -1.88
4 fixed <NA> tratamentoStrips 25% -0.773 1.04 -0.742
5 fixed <NA> tratamentoStrips 50% -1.15 1.04 -1.1
6 fixed <NA> ano -0.346 0.121 -2.86
7 fixed <NA> tratamentoNuclei 50%:ano 0.022 0.171 0.13
8 fixed <NA> tratamentoPlanting 100%:ano 0.136 0.171 0.794
9 fixed <NA> tratamentoStrips 25%:ano 0.008 0.171 0.045
10 fixed <NA> tratamentoStrips 50%:ano 0.067 0.171 0.389
11 ran_pars bloco sd__(Intercept) 0.298 NA NA
12 ran_pars Residual sd__Observation 0.383 NA NA
fig_data_cost_raw <- cost_raw |>
dplyr::group_by(ano, tratamento)|>
dplyr::summarise(
cost_cc = mean(cc_cost),
SE_cc = sd(cc_cost)/sqrt(dplyr::n()),
cost_agb = mean(agb_cost),
SE_agb = sd(agb_cost)/sqrt(dplyr::n()),
.groups="drop")|>
dplyr::mutate(tratamento = ordered(tratamento, levels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%")),
tratamento = forcats::fct_recode(tratamento,
"Nuclei\n25%" = "Nuclei 25%",
"Strips\n25%" = "Strips 25%",
"Nuclei\n50%" = "Nuclei 50%",
"Strips\n50%" = "Strips 50%",
"Planting\n100%" = "Planting 100%"
))letras <- data.frame(
tratamento = rep(levels(as.factor(fig_data_cost_raw$tratamento)), each=2),
ano = c(5, 7),
label = c("B "," a",
"AB "," a",
"AB "," a",
"AB "," a",
"A "," a"),
y = c(33, 16,
18, 10,
18, 8,
18, 10,
10, 6))
fig_cost1 <- fig_data_cost_raw |>
ggplot(aes(x = tratamento, y = cost_cc, fill = as.factor(ano))) +
geom_bar(position = position_dodge2(preserve = "single"),
stat = "identity", color = "black") +
geom_errorbar(aes(ymin = cost_cc - SE_cc, ymax = cost_cc + SE_cc), width = .2, position = position_dodge(0.9)) +
scale_fill_manual(
labels = c('5 years', '7 years'),
values = c("#CAF0F8", "#00B4D8")) +
geom_text(data = letras, aes(x = tratamento, y = y, label = label), size = 3.5, color="black", vjust = -0.5) +
scale_y_continuous(expand=c(0,0), limits=c(0,38), breaks = c(0, 10,20, 30))+
labs(fill = "",
y = expression("CE Canopy cover (" * "\u0024" * ".%"^{-1} * ")"),
x = "") +
theme_classic()+
theme(
axis.text.x = element_blank(),
axis.text.y = element_text(color="black", size=10),
axis.title.y = element_text(color="black", size=10),
plot.margin = margin(t=5, 0, 0, 0, "mm"),
legend.position = c(0.85, 0.95))
fig_cost1letras <- data.frame(
tratamento = rep(levels(as.factor(fig_data_cost_raw$tratamento)), each=2),
ano = c(5, 7),
label = c("B "," b",
"A "," a",
"AB "," ab",
"A "," ab",
"A "," a"),
y = c(210, 100, 53, 25, 53, 30, 53, 33, 30, 20)
)
fig_cost2 <- fig_data_cost_raw |>
ggplot(aes(x = tratamento, y = cost_agb, fill = as.factor(ano))) +
geom_bar(position = position_dodge2(preserve = "single"),
stat = "identity", color = "black") +
geom_errorbar(aes(ymin = cost_agb - SE_agb, ymax = cost_agb + SE_agb),
width = .2, position = position_dodge(0.9)) +
scale_fill_manual(
labels = c('5 years', '7 years'),
values = c("#CAF0F8", "#00B4D8")) +
geom_text(data = letras, aes(x = tratamento, y = y, label = label), size = 3.5, color="black", vjust = -0.5) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 240), breaks = c(0,50,100, 150, 200))+
labs(fill = "",
y = expression("CE AGB (" * "\u0024" * ".Mg"^{-1} * ")"),
x = "") +
theme_classic()+
theme(
axis.line.x = element_line(color = "black", linewidth = 0.8, linetype = "solid"),
axis.line.y = element_line(color = "black", linewidth = 0.8, linetype = "solid"),
axis.ticks = element_line(color = "black", linewidth = 0.5),
axis.text.x = element_text(color="black", hjust = .5, size=10),
axis.text.y = element_text(color="black", size=10),
legend.position = 'none',
axis.title.y = element_text(color="black", size=10))
fig_cost2fig_cost <- cowplot::plot_grid(fig_cost1, fig_cost2, ncol = 1,
align = "v", axis = "l",
labels = c("(a)", "(b)"), label_size = 8,
label_x = 0.05, label_y = c(1,1.1))
fig_cost# ggsave("fig_cost.png", plot=fig_cost, width=90, height = 120 ,unit="mm")fig_data_cost_raw |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 2))) # A tibble: 10 × 6
ano tratamento cost_cc SE_cc cost_agb SE_agb
<dbl> <ord> <dbl> <dbl> <dbl> <dbl>
1 5 "Nuclei\n25%" 23.3 9.43 133. 69.9
2 5 "Nuclei\n50%" 15.5 1.08 46.4 4.84
3 5 "Planting\n100%" 7.07 1.09 25.6 3.28
4 5 "Strips\n25%" 12.8 1.8 44.6 5.95
5 5 "Strips\n50%" 14.5 3.44 43.3 8.74
6 7 "Nuclei\n25%" 10.5 4.82 58.6 26.0
7 7 "Nuclei\n50%" 6.26 0.94 24.4 2.88
8 7 "Planting\n100%" 5.35 0.11 16.6 1.47
9 7 "Strips\n25%" 7.52 1.37 22.6 2.66
10 7 "Strips\n50%" 6.63 1.92 26.3 7.65